#include <math.h>
#include "floatimage.h"
#include "stdio.h"
#include "memory.h"
#include <malloc.h>
#include "bmp.h"
#include "gamma.h"
#include "rad_color.h"
#include "wvtimage.h"

#ifdef FLOATIMAGE_RGBEFILES
	#include "rad_color.h"
	#include "rad_header.h"
	#include "rad_resolu.h"
#endif

// floating point TIFF format headers
#ifdef FLOATIMAGE_USETIFF
#include "tiffio.h"
#endif

#ifdef FLOATIMAGE_USEJPEG
#include "jpegfile.h"
#endif

#ifndef UINT
#define UINT unsigned 
#endif
#ifndef BYTE
#define BYTE char
#endif
#ifndef TRUE
#define TRUE 1
#define FALSE 0
#endif

#define REVERSE2BYTES(source, dest)		\
{										\
	char *j = (char *) source;			\
	char *dj = (char *) dest;			\
	dj[0] = j[1];						\
	dj[1] = j[0];						\
}


#define REVERSEBYTES(source, dest)		\
{										\
	char *j = (char *) source;			\
	char *dj = (char *) dest;			\
	dj[0] = j[3];						\
	dj[1] = j[2];						\
	dj[2] = j[1];						\
	dj[3] = j[0];						\
}

#define REVERSE2BYTES_INPLACE(sourcedest)		\
{												\
	char *j = (char *) sourcedest;				\
	j[0] ^= j[1]; j[1] ^= j[0]; j[0] ^= j[1];	\
}

#define REVERSEBYTES_INPLACE(sourcedest)		\
{												\
	char *j = (char *) sourcedest;				\
	j[0] ^= j[3]; j[3] ^= j[0]; j[0] ^= j[3];	\
	j[1] ^= j[2]; j[2] ^= j[1]; j[1] ^= j[2];	\
}


#define  subcolor(c1,c2)	((c1)[0]-=(c2)[0],(c1)[1]-=(c2)[1],(c1)[2]-=(c2)[2])
#define  lumcolor(c1)	(0.3f*(c1)[0]+0.6f*(c1)[1]+0.1f*(c1)[2])


#ifndef M_PI
# define M_PI	3.14159265358979323846264338327950288
#endif

#define m_boxvariance(n) (((2*(n)*(n)*(n)+3*(n)*(n)+(n))/3.f)*(4*(n)+2)/((2*(n)+1)*(2.f*(n)+1)))

inline unsigned short int float2short(float f)
{
  unsigned int bits = * (unsigned int *) &f;
  return (bits >> 15);
}


inline unsigned short int signedfloat2short(float f)
{
  unsigned int bits = * (unsigned int *) &f;
  return (bits >> 16);
}


inline float short2float(unsigned short int s)
{
	unsigned int bits = s << 15;
	return * (float *) &bits;
}

inline float short2signedfloat(unsigned short int s)
{
	unsigned int bits = s << 16;
	return * (float *) &bits;
}


void floatimage::ClipLower(float min, float replace)
{
	float *f = &line(0)[0].r;
	float *last = &line(sizey-1)[sizex-1].b;

	while (f <= last)
	{
		float r = f[0];
		float g = f[1];
		float b = f[2];

		if ((f[0] < min) && (f[1] < min) && (f[2] < min))
		{
			f[0] = replace;
			f[1] = replace;
			f[2] = replace;
		}

		f += 3;
	}
}

#ifdef FLOATIMAGE_RGBEFILES

int floatimage::LoadHDR(const char *filename)
{
    FILE *infile = fopen(filename, "rb");
	int new_sizex,new_sizey;

	sprintf(name, "%s", filename);

	if (infile == NULL)
		return 0;

	getheader(infile);

	int mode = fgetresolu( &new_sizex, &new_sizey, infile);

	Allocate(new_sizex, new_sizey);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	for (int y = 0; y < sizey; y++)
	{
		if (freadscan((COLOR *) line(y), sizex, infile) < 0) 
		{ 
			Release();
			fclose(infile);
			return 0;
		}
	}

	fclose(infile);

	if (mode != (YMAJOR | YDECR))
		return 2;

	return 1;
}

int floatimage::SaveHDR(const char *filename)
{
	FILE *outfile = fopen(filename, "wb");

    fprintf(outfile, "#?RADIANCE");
	fputc(10, outfile);
    fprintf(outfile, "# %s", "Made with 100% pure HDR Shop");
	fputc(10, outfile);
    fprintf(outfile, "FORMAT=32-bit_rle_rgbe");
	fputc(10, outfile);
    fprintf(outfile, "EXPOSURE=%25.13f", 1.0);
	fputc(10, outfile);
	fputc(10, outfile);
    fprintf(outfile, "-Y %d +X %d", sizey, sizex);
	fputc(10, outfile);
	
	for (int y = 0; y < sizey; y++)
	{
		if (fwritescan((COLOR *) &line(y)[0].r, sizex, outfile) < 0)
		{
			// error
			fclose(outfile);
			return 0;
		}
	}

	fclose(outfile);

	return 1;
}

#endif

void floatimage::SwapByteOrder()
{
	for (int y = 0; y < sizey; y++)
		for (int x = 0; x < sizex; x++)
		{
			REVERSEBYTES_INPLACE(&line(y)[x].r);
			REVERSEBYTES_INPLACE(&line(y)[x].g);
			REVERSEBYTES_INPLACE(&line(y)[x].b);
		}
}


int floatimage::LoadPFM(const char *filename)
{
	FILE *infile = fopen(filename, "rb");

	sprintf(name, "%s", filename);

	unsigned int a = fgetc(infile);
	int b = fgetc(infile);
	int junk = fgetc(infile);

	if ((a != 'P') || ((b != 'F') && (b != 'f')))
	{
		fclose(infile);
		return 0;
	}

	b = (b == 'F');		// 'F' = RGB,  'f' = monochrome

	int x, y;
	fscanf(infile, "%d %d%c", &x, &y, &junk);
	if ((x <= 0) || (y <= 0))
	{
		fclose(infile);
		return 0;
	}

	float scalefactor;
	fscanf(infile, "%f%c", &scalefactor, &junk);


	Allocate(x, y);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	a = sizex * (b ? 3 : 1);
	float *fbuf = new float[a];
	fRGB *cur = &line(0)[0];
	for (y = 0; y < sizey; y++)
	{
		if (fread(fbuf, sizeof(float), a, infile) != a)
		{
			fclose(infile);
			delete fbuf;
			return 0;
		}

		float *temp = fbuf;
		for (x = 0; x < sizex; x++)
		{
			if (b)			// color
			{
				if (scalefactor > 0.0)		// MSB
				{
					REVERSEBYTES(temp++, &cur->r);
					REVERSEBYTES(temp++, &cur->g);
					REVERSEBYTES(temp++, &cur->b);
				}
				else						// LSB
				{
					cur->r = *temp++;
					cur->g = *temp++;
					cur->b = *temp++;
				}
			}
			else			// black and white
			{
				float c;
				if (scalefactor > 0.0)		// MSB
				{
					REVERSEBYTES(temp++, &c);
				}
				else						// LSB
					c = *temp++;
				cur->r = cur->g = cur->b = c;
			}
			cur++;
		}
	}

	delete fbuf;
	fclose(infile);

	FlipVertical();

	return 1;
}

int floatimage::LoadPFMCropped(const char *filename, int startX, int startY, int width, int height)
{
	FILE *infile = fopen(filename, "rb");

	if (!infile)
		return 0 ;

	sprintf(name, "%s", filename);

	unsigned int a = fgetc(infile);
	int b = fgetc(infile);
	int junk = fgetc(infile);

	if ((a != 'P') || ((b != 'F') && (b != 'f')))
	{
		fclose(infile);
		return 0;
	}

	b = (b == 'F');		// 'F' = RGB,  'f' = monochrome

	int fullX, fullY;
	fscanf(infile, "%d %d%c", &fullX, &fullY, &junk);

	// Invalid PFM image size...
	if ((fullX <= 0) || (fullY <= 0))
	{
		fclose(infile);
		return 0;
	}

	float scalefactor;
	fscanf(infile, "%f%c", &scalefactor, &junk);
	

	// First, do all the checks to insure we're safe...
	if  ((startX < 0) ||
		 (startY < 0) ||
		 ((fullX-startX) <= 0) || 
		 ((fullY-startY) <= 0))
	{
		fclose(infile) ;
		return 0 ;
	}

	// Second, figure out where the third '0x0A' is 'cuz that's where you need to read from...
	fseek(infile, 0, SEEK_SET) ;

	int oaCount = 0 ;
	int startOffset = 0 ;
	while(oaCount < 3)
	{
		unsigned char bob = fgetc(infile) ;
		if (bob == 0x0A)
			oaCount++ ;

		startOffset++ ;
	}

	// Looks good to me, load the file now yo...
	Allocate(width, height) ;

	// Can't allocate enough memory?
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	a = width * (b ? 3 : 1);

	float *fbuf = new float[a];
	fRGB *cur = &line(0)[0];

	int x, y ;

	// Reads it in upside-down so you don't need the flipVertical at the end...
	for (y = fullY-startY-1; y >= fullY-height-startY; y--)
	{
		fseek(infile, (((y*fullX*(b ? 3 : 1)) + (startX*(b ? 3 : 1))) * sizeof(float)) + startOffset, SEEK_SET) ;

		if (fread(fbuf, sizeof(float), a, infile) != a)
		{
			// Something's amiss with fread
			fclose(infile);
			delete fbuf;
			return 0;
		}

		float *temp = fbuf;
		for (x = 0; x < sizex; x++)
		{
			if (b)			// color
			{
				if (scalefactor > 0.0)		// MSB
				{
					REVERSEBYTES(temp++, &cur->r);
					REVERSEBYTES(temp++, &cur->g);
					REVERSEBYTES(temp++, &cur->b);
				}
				else						// LSB
				{
					cur->r = *temp++;
					cur->g = *temp++;
					cur->b = *temp++;
				}
			}
			else			// black and white
			{
				float c;
				if (scalefactor > 0.0)		// MSB
				{
					REVERSEBYTES(temp++, &c);
				}
				else						// LSB
					c = *temp++;
				cur->r = cur->g = cur->b = c;
			}
			cur++;
		}
	}

	delete fbuf;
	fclose(infile);

//	FlipVertical();

	return 1 ;
}


int floatimage::SavePFM(const char *filename)
{
	FILE *outfile = fopen(filename, "wb");

	if (!outfile)
		return 0;

	fputc('P', outfile);
	fputc('F', outfile);
	fputc(0x0a, outfile);

	fprintf(outfile, "%d %d", sizex, sizey);
	fputc(0x0a, outfile);

	fprintf(outfile, "%f", -1.0f);
	fputc(0x0a, outfile);

	int result = 1;
	for (int y = sizey - 1; y >= 0; y--)
	{
		if (fwrite(&line(y)[0].r, sizeof(float) * sizex, 3, outfile) < 3)
			result = 0;
	}

	fclose(outfile);

	return result;
}


int floatimage::SavePFM_channel(const char *filename, int channel)
{
	FILE *outfile = fopen(filename, "wb");

	if (!outfile)
		return 0;

	fputc('P', outfile);
	fputc('f', outfile);
	fputc(0x0a, outfile);

	fprintf(outfile, "%d %d", sizex, sizey);
	fputc(0x0a, outfile);

	fprintf(outfile, "%f", -1.0f);
	fputc(0x0a, outfile);


	int result = 1;
	for (int y = sizey - 1; y >= 0; y--)
	{
		float *cur = line(y)[0].channel(channel);
		for (int x = 0; x < sizex; x++)
		{
			if (fwrite(cur, sizeof(float), 1, outfile) < 1)
				result = 0;
			cur += 3;
		}
	}

	fclose(outfile);

	return result;
}



int floatimage::SaveRAW(const char *filename)
{
	FILE *outfile = fopen(filename, "wb");
	if (!outfile)
	{
		return 0;
	}

	for (int y = sizey - 1; y >= 0; y--)
	{
		fwrite(&line(y)[0].r, sizeof(float), 3 * sizex, outfile);
	}

	fclose(outfile);
	return 1;
}

int floatimage::LoadRAWBytes(const char *filename, int width, int height, int rgb, float gamma, int header)
{
	InitGamma(gamma);
	
	FILE *infile = fopen(filename, "rb");
	if (!infile)
	{
		return 0;
	}

	fseek(infile, header, SEEK_CUR);

	Allocate(width, height);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	unsigned int numbytes = (rgb ? 3 * width : width);
	unsigned char *fbuf = new unsigned char[numbytes];

	fRGB *cur = &line(0)[0];
	for (int y = 0; y < height; y++)
	{
		if (fread(fbuf, sizeof(unsigned char), numbytes, infile) != numbytes)
		{
			delete fbuf;
			fclose(infile);
//			MessageBox(NULL, "Error Reading Image File", blah->name, MB_OK);
			return 0;
		}

		unsigned char *temp = fbuf;
		for (int x = 0; x < width; x++)
		{
			if (rgb)
			{
				cur->r = gamma_table[*temp++];
				cur->g = gamma_table[*temp++];
				cur->b = gamma_table[*temp++];
			}
			else
			{
				cur->r = cur->g = cur->b = gamma_table[*temp++];
			}

			cur ++;
		}
	}

	delete fbuf;
	fclose(infile);

	FlipVertical();

	sprintf(name, "%s", filename);

	return 1;
}

int floatimage::LoadRAWShorts(const char *filename, int width, int height, int rgb, int byte_reverse, int header)
{
	FILE *infile = fopen(filename, "rb");
	if (!infile)
	{
		return 0;
	}

	fseek(infile, header, SEEK_CUR);
	
	Allocate(width, height);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	unsigned int numshorts = (rgb ? 3 * width : width);
	short int *fbuf = new short int[numshorts];

	fRGB *cur = &line(0)[0];
	for (int y = 0; y < height; y++)
	{
		if (fread(fbuf, sizeof(short int), numshorts, infile) != numshorts)
		{
			delete fbuf;
			fclose(infile);
//			MessageBox(NULL, "Error Reading Image File", blah->name, MB_OK);
			return 0;
		}

		short int *temp = fbuf;
		for (int x = 0; x < width; x++)
		{
			if (rgb)
			{
				if (byte_reverse)		// MSB
				{
					REVERSE2BYTES(temp++, &cur->r);
					REVERSE2BYTES(temp++, &cur->g);
					REVERSE2BYTES(temp++, &cur->b);
				}
				else
				{
					cur->r = *temp++;
					cur->g = *temp++;
					cur->b = *temp++;
				}
			}
			else
			{
				short c;
				if (byte_reverse)		// MSB
				{
					REVERSE2BYTES(temp++, &c);
				}
				else						// LSB
					c = *temp++;
				cur->r = cur->g = cur->b = c;
			}

			cur ++;
		}
	}

	delete fbuf;
	fclose(infile);

	FlipVertical();

	sprintf(name, "%s", filename);

	return 1;
}


int floatimage::LoadRAWUnsignedShorts(const char *filename, int width, int height, int rgb, int byte_reverse, int header)
{
	FILE *infile = fopen(filename, "rb");
	if (!infile)
	{
		return 0;
	}
	
	fseek(infile, header, SEEK_CUR);

	Allocate(width, height);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	unsigned int numshorts = (rgb ? 3 * width : width);
	unsigned short int *fbuf = new unsigned short int[numshorts];

	fRGB *cur = &line(0)[0];
	for (int y = 0; y < height; y++)
	{
		if (fread(fbuf, sizeof(unsigned short int), numshorts, infile) != numshorts)
		{
			delete fbuf;
			fclose(infile);
			return 0;
		}

		unsigned short int *temp = fbuf;
		for (int x = 0; x < width; x++)
		{
			if (rgb)
			{
				if (byte_reverse)		// MSB
				{
					REVERSE2BYTES(temp++, &cur->r);
					REVERSE2BYTES(temp++, &cur->g);
					REVERSE2BYTES(temp++, &cur->b);
				}
				else
				{
					cur->r = *temp++;
					cur->g = *temp++;
					cur->b = *temp++;
				}
			}
			else
			{
				unsigned short c;
				if (byte_reverse)		// MSB
				{
					REVERSE2BYTES(temp++, &c);
				}
				else						// LSB
					c = *temp++;
				cur->r = cur->g = cur->b = c;
			}

			cur ++;
		}
	}

	delete fbuf;
	fclose(infile);

	FlipVertical();

	sprintf(name, "%s", filename);

	return 1;
}


int floatimage::LoadRAWLongs(const char *filename, int width, int height, int rgb, int byte_reverse, int header)
{
	FILE *infile = fopen(filename, "rb");
	if (!infile)
	{
		return 0;
	}

	fseek(infile, header, SEEK_CUR);
	
	Allocate(width, height);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	unsigned int numlongs = (rgb ? 3 * width : width);
	int *fbuf = new int[numlongs];

	fRGB *cur = &line(0)[0];
	for (int y = 0; y < height; y++)
	{
		if (fread(fbuf, sizeof(int), numlongs, infile) != numlongs)
		{
			delete fbuf;
			fclose(infile);
			return 0;
		}

		int *temp = fbuf;
		for (int x = 0; x < width; x++)
		{
			if (rgb)
			{
				if (byte_reverse)		// MSB
				{
					REVERSEBYTES(temp++, &cur->r);
					REVERSEBYTES(temp++, &cur->g);
					REVERSEBYTES(temp++, &cur->b);
				}
				else
				{
					cur->r = (float)*temp++;
					cur->g = (float)*temp++;
					cur->b = (float)*temp++;
				}
			}
			else
			{
				unsigned long c;
				if (byte_reverse)		// MSB
				{
					REVERSEBYTES(temp++, &c);
				}
				else						// LSB
					c = *temp++;
				cur->r = cur->g = cur->b = (float)c;
			}

			cur ++;
		}
	}

	delete fbuf;
	fclose(infile);

	FlipVertical();

	sprintf(name, "%s", filename);

	return 1;
}


int floatimage::LoadRAWUnsignedLongs(const char *filename, int width, int height, int rgb, int byte_reverse, int header)
{
	FILE *infile = fopen(filename, "rb");
	if (!infile)
	{
		return 0;
	}

	fseek(infile, header, SEEK_CUR);
	
	Allocate(width, height);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	unsigned int numlongs = (rgb ? 3 * width : width);
	unsigned int *fbuf = new unsigned int[numlongs];

	fRGB *cur = &line(0)[0];
	for (int y = 0; y < height; y++)
	{
		if (fread(fbuf, sizeof(unsigned int), numlongs, infile) != numlongs)
		{
			delete fbuf;
			fclose(infile);
			return 0;
		}

		unsigned int *temp = fbuf;
		for (int x = 0; x < width; x++)
		{
			if (rgb)
			{
				if (byte_reverse)		// MSB
				{
					REVERSEBYTES(temp++, &cur->r);
					REVERSEBYTES(temp++, &cur->g);
					REVERSEBYTES(temp++, &cur->b);
				}
				else
				{
					cur->r = (float)*temp++;
					cur->g = (float)*temp++;
					cur->b = (float)*temp++;
				}
			}
			else
			{
				unsigned long c;
				if (byte_reverse)		// MSB
				{
					REVERSEBYTES(temp++, &c);
				}
				else						// LSB
					c = *temp++;
				cur->r = cur->g = cur->b = (float)c;
			}

			cur ++;
		}
	}

	delete fbuf;
	fclose(infile);

	FlipVertical();

	sprintf(name, "%s", filename);

	return 1;
}


int floatimage::LoadRAW(const char *filename, int width, int height, int rgb, int byte_reverse, int header)
{
	FILE *infile = fopen(filename, "rb");
	if (!infile)
	{
		return 0;
	}

	fseek(infile, header, SEEK_CUR);
	
	Allocate(width, height);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	unsigned int numfloats = (rgb ? 3 * width : width);
	float *fbuf = new float[numfloats];

	fRGB *cur = &line(0)[0];
	for (int y = 0; y < height; y++)
	{
		if (fread(fbuf, sizeof(float), numfloats, infile) != numfloats)
		{
			delete fbuf;
			fclose(infile);
			return 0;
		}

		float *temp = fbuf;
		for (int x = 0; x < width; x++)
		{
			if (rgb)
			{
				if (byte_reverse)		// MSB
				{
					REVERSEBYTES(temp++, &cur->r);
					REVERSEBYTES(temp++, &cur->g);
					REVERSEBYTES(temp++, &cur->b);
				}
				else						// LSB
				{
					cur->r = *temp++;
					cur->g = *temp++;
					cur->b = *temp++;
				}
			}
			else
			{
				float c;
				if (byte_reverse)		// MSB
				{
					REVERSEBYTES(temp++, &c);
				}
				else						// LSB
					c = *temp++;
				cur->r = cur->g = cur->b = c;
			}

			cur ++;
		}
	}

	delete fbuf;
	fclose(infile);

	sprintf(name, "%s", filename);

	return 1;
}

#ifdef FLOATIMAGE_USETIFF
int floatimage::LoadTIFSize(const char *filename)
{
	long ImageWidth = 0;
	long ImageLength = 0;

	TIFF *tiffImage = TIFFOpen(filename, "ru");

	if (!tiffImage) 
		return 0;

	TIFFGetField(tiffImage, TIFFTAG_IMAGEWIDTH, &ImageWidth);
	TIFFGetField(tiffImage, TIFFTAG_IMAGELENGTH, &ImageLength);

	TIFFClose(tiffImage);

	sizex = ImageWidth;
	sizey = ImageLength;

	return 1;
}



#ifdef  NTSC
#define  CIE_x_r		0.670f		// standard NTSC primaries
#define  CIE_y_r		0.330f
#define  CIE_x_g		0.210f
#define  CIE_y_g		0.710f
#define  CIE_x_b		0.140f
#define  CIE_y_b		0.080f
#define  CIE_x_w		0.3333f		// use true white
#define  CIE_y_w		0.3333f
#else
#define  CIE_x_r		0.640f		// nominal CRT primaries
#define  CIE_y_r		0.330f
#define  CIE_x_g		0.290f
#define  CIE_y_g		0.600f
#define  CIE_x_b		0.150f
#define  CIE_y_b		0.060f
#define  CIE_x_w		0.3333f		// use true white
#define  CIE_y_w		0.3333f
#endif

COLORMAT  rgb2xyzmat =				// RGB to XYZ
{
	{CIE_x_r * CIE_C_rD / CIE_D, CIE_x_g * CIE_C_gD / CIE_D, CIE_x_b * CIE_C_bD / CIE_D},
	{CIE_y_r * CIE_C_rD / CIE_D, CIE_y_g * CIE_C_gD / CIE_D, CIE_y_b * CIE_C_bD / CIE_D},
	{(1.f - CIE_x_r - CIE_y_r) * CIE_C_rD / CIE_D,
	 (1.f - CIE_x_g - CIE_y_g) * CIE_C_gD / CIE_D,
	 (1.f - CIE_x_b - CIE_y_b) * CIE_C_bD / CIE_D}
};

COLORMAT  xyz2rgbmat =				/* XYZ to RGB */
{
	{(CIE_y_g - CIE_y_b - CIE_x_b*CIE_y_g + CIE_y_b*CIE_x_g)/CIE_C_rD,
	 (CIE_x_b - CIE_x_g - CIE_x_b*CIE_y_g + CIE_x_g*CIE_y_b)/CIE_C_rD,
	 (CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g)/CIE_C_rD},
	{(CIE_y_b - CIE_y_r - CIE_y_b*CIE_x_r + CIE_y_r*CIE_x_b)/CIE_C_gD,
	 (CIE_x_r - CIE_x_b - CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r)/CIE_C_gD,
	 (CIE_x_b*CIE_y_r - CIE_x_r*CIE_y_b)/CIE_C_gD},
	{(CIE_y_r - CIE_y_g - CIE_y_r*CIE_x_g + CIE_y_g*CIE_x_r)/CIE_C_bD,
	 (CIE_x_g - CIE_x_r - CIE_x_g*CIE_y_r + CIE_x_r*CIE_y_g)/CIE_C_bD,
	 (CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r)/CIE_C_bD}
};


int floatimage::LoadTIF_NoError(const char *filename, float image_gamma, int * bits)
{

	int result = LoadTIF(filename, image_gamma, bits);
	return result;
}


int floatimage::LoadTIF(const char *filename, float image_gamma, int * bits)
{
	long ImageWidth = 0;
	long ImageLength = 0;
	long BitsPerSample = 0;
	long SamplesPerPixel = 0;

	int CHANNELS;
	int PIXELS;

	int i, j;
	TIFF *tiffImage;

	tiffImage = TIFFOpen(filename, "ru");

	if (!tiffImage) 
		return 0;

	TIFFGetField(tiffImage, TIFFTAG_IMAGEWIDTH, &ImageWidth);
	TIFFGetField(tiffImage, TIFFTAG_IMAGELENGTH, &ImageLength);
	TIFFGetField(tiffImage, TIFFTAG_SAMPLESPERPIXEL, &SamplesPerPixel);
	TIFFGetField(tiffImage, TIFFTAG_BITSPERSAMPLE, &BitsPerSample);

	long comp = 0;

	if (!TIFFGetField(tiffImage, TIFFTAG_COMPRESSION, &comp))
        comp = COMPRESSION_NONE;

	long phot = 0;

	int convertXYZ = 0;

	if (TIFFGetFieldDefaulted(tiffImage, TIFFTAG_PHOTOMETRIC, &phot))
	{
		switch (phot)
		{
			case PHOTOMETRIC_LOGLUV:
			case PHOTOMETRIC_LOGL:   /* include only if prepared for Y-only data */
					if ((comp != COMPRESSION_SGILOG) && (comp != COMPRESSION_SGILOG24))
							return 0;	// quiterr("Only support SGILOG compressed LogLuv data");
					TIFFSetField(tiffImage, TIFFTAG_SGILOGDATAFMT, SGILOGDATAFMT_FLOAT);
					BitsPerSample = 32;
					convertXYZ = 1;
					break;
			/* ...other supported photometric types */
			default:
//					return 0; // quiterr("TIFF has unsupported photometric type");
				break;
		}
	}

	CHANNELS = SamplesPerPixel;
	PIXELS = ImageWidth * ImageLength;

//	if ((CHANNELS != 1) && (CHANNELS != 3))
//		return 0;

	if ((BitsPerSample != 32) && (BitsPerSample != 8) && (BitsPerSample != 16))
		return 0;

	Allocate(ImageWidth, ImageLength);

	if (!buf)
		return 0;

	if (bits)
		*bits = BitsPerSample;
    
	if (BitsPerSample == 32)		// floats
	{
		float *scanline = (float *) _TIFFmalloc(TIFFScanlineSize(tiffImage));
		_TIFFmemset(scanline, 0, TIFFScanlineSize(tiffImage));

		for (i = 0; i < sizey; i++)
		{
//			TIFF_Error_Flag = 0;
			TIFFReadScanline(tiffImage, scanline, i, 0);
//			if (TIFF_Error_Flag)
//				break;

			float *buf = scanline;
			fRGB *linei = line(i);

			if (CHANNELS==1)
			{
				for (j = 0; j < sizex; j++)
				{
					float c = *buf++;
					linei->Set(c, c, c);
					linei++;
				}
			}
			else
			{
				for (j = 0; j < sizex; j++)
				{
					linei->Set(buf[0], buf[1], buf[2]);
					linei++;
					buf += CHANNELS;
				} 
			}
		}
		_TIFFfree((void *) scanline);

		if (convertXYZ)
		{
			// convert
		}
	}
	else if (BitsPerSample == 16)	// 16 bit linear
	{
		unsigned char *scanline = (unsigned char *) _TIFFmalloc(TIFFScanlineSize(tiffImage));
		_TIFFmemset(scanline, 0, TIFFScanlineSize(tiffImage));

		InitGamma(image_gamma);

		for (i = 0; i < sizey; i++)
		{
//			TIFF_Error_Flag = 0;
			TIFFReadScanline(tiffImage, scanline, i, 0);
//			if (TIFF_Error_Flag)
//				break;

			unsigned short int *buf = (unsigned short int *) scanline;
			fRGB *linei = line(i);

			if (CHANNELS==1)
			{
				for (j = 0; j < sizex; j++)
				{
					float c = (*buf++) / 65536.0f;
					linei->Set(c, c, c);
					linei++;
				}
			}
			else
			{
				for (j = 0; j < sizex; j++)
				{
					linei->Set(buf[0] / 65536.0f, buf[1] / 65536.0f, buf[2] / 65536.0f);
					linei++;
					buf += CHANNELS;
				} 
			}
		}
		_TIFFfree((void *) scanline);
	}  
	else if (BitsPerSample == 8)	// bytes
	{
		unsigned char *scanline = (unsigned char *) _TIFFmalloc(TIFFScanlineSize(tiffImage));
		_TIFFmemset(scanline, 0, TIFFScanlineSize(tiffImage));

		InitGamma(image_gamma);

		for (i = 0; i < sizey; i++)
		{
//			TIFF_Error_Flag = 0;
			TIFFReadScanline(tiffImage, scanline, i, 0);
//			if (TIFF_Error_Flag)
//				break;

			unsigned char *buf = scanline;
			fRGB *linei = line(i);

			if (CHANNELS==1)
			{
				for (j = 0; j < sizex; j++)
				{
					float c = gamma_table[*buf++];
					linei->Set(c, c, c);
					linei++;
				}
			}
			else
			{
				for (j = 0; j < sizex; j++)
				{
					linei->Set(gamma_table[buf[0]], gamma_table[buf[1]], gamma_table[buf[2]]);
					linei++;
					buf += CHANNELS;
				} 
			}
		}
		_TIFFfree((void *) scanline);
	}  

	TIFFClose(tiffImage);

	if (convertXYZ)
	{
		for (int y = 0; y < sizey; y++)
			for (int x = 0; x < sizex; x++)
			{
				double max = line(y)[x].g;
	// TODO FIX THIS BUGGER			line(y)[x].transform(xyz2rgbmat);
//				line(y)[x].clip_gamut_lower(max);
			}
	}

	sprintf(name, "%s", filename);

	return 1;
}


int floatimage::LoadTIF_Curve(const char *filename, float rcurve[256], float gcurve[256], float bcurve[256])
{
	long ImageWidth = 0;
	long ImageLength = 0;
	long BitsPerSample = 0;
	long SamplesPerPixel = 0;

	int CHANNELS;
	int PIXELS;

	int i, j;
	TIFF *tiffImage;

	tiffImage = TIFFOpen(filename, "ru");

	if (!tiffImage) 
		return 0;

	TIFFGetField(tiffImage, TIFFTAG_IMAGEWIDTH, &ImageWidth);
	TIFFGetField(tiffImage, TIFFTAG_IMAGELENGTH, &ImageLength);
	TIFFGetField(tiffImage, TIFFTAG_SAMPLESPERPIXEL, &SamplesPerPixel);
	TIFFGetField(tiffImage, TIFFTAG_BITSPERSAMPLE, &BitsPerSample);

	long comp = 0;

	if (!TIFFGetField(tiffImage, TIFFTAG_COMPRESSION, &comp))
        comp = COMPRESSION_NONE;

	long phot = 0;

	int convertXYZ = 0;

	if (TIFFGetFieldDefaulted(tiffImage, TIFFTAG_PHOTOMETRIC, &phot))
	{
		convertXYZ = 1;
	//	switch (phot)
	//	{
//			case PHOTOMETRIC_LOGLUV:
//			case PHOTOMETRIC_LOGL:   /* include only if prepared for Y-only data */
//					if (comp != COMPRESSION_SGILOG && comp != COMPRESSION_SGILOG24)
//							return 0;	// quiterr("Only support SGILOG compressed LogLuv data");
//					TIFFSetField(tiffImage, TIFFTAG_SGILOGDATAFMT, SGILOGDATAFMT_FLOAT);
	//				break;
			/* ...other supported photometric types */
	//		default:
//					return 0; // quiterr("TIFF has unsupported photometric type");
	//			break;
	//	}
	}

	CHANNELS = SamplesPerPixel;
	PIXELS = ImageWidth * ImageLength;

//	if ((CHANNELS != 1) && (CHANNELS != 3))
//		return 0;

	if ((BitsPerSample != 32) && (BitsPerSample != 8) && (BitsPerSample != 16))
		return 0;

	Allocate(ImageWidth, ImageLength);

	if (!buf)
		return 0;
    
	if (BitsPerSample == 32)		// floats
	{
		float *scanline = (float *) _TIFFmalloc(TIFFScanlineSize(tiffImage));
		_TIFFmemset(scanline, 0, TIFFScanlineSize(tiffImage));

		for (i = 0; i < sizey; i++)
		{
//			TIFF_Error_Flag = 0;
			TIFFReadScanline(tiffImage, scanline, i, 0);
//			if (TIFF_Error_Flag)
//				break;

			float *buf = scanline;
			fRGB *linei = line(i);

			if (CHANNELS==1)
			{
				for (j = 0; j < sizex; j++)
				{
					float c = *buf++;
					linei->Set(c, c, c);
					linei++;
				}
			}
			else
			{
				for (j = 0; j < sizex; j++)
				{
					linei->Set(buf[0], buf[1], buf[2]);
					linei++;
					buf += CHANNELS;
				} 
			}
		}
		_TIFFfree((void *) scanline);

		if (convertXYZ)
		{
			// convert
		}
	}
	else if (BitsPerSample == 16)	// 16 bit curve
	{
		unsigned char *scanline = (unsigned char *) _TIFFmalloc(TIFFScanlineSize(tiffImage));
		_TIFFmemset(scanline, 0, TIFFScanlineSize(tiffImage));

		for (i = 0; i < sizey; i++)
		{
//			TIFF_Error_Flag = 0;
			TIFFReadScanline(tiffImage, scanline, i, 0);
//			if (TIFF_Error_Flag)
//				break;

			unsigned short int *buf = (unsigned short int *) scanline;
			fRGB *linei = line(i);

			if (CHANNELS==1)
			{
				for (j = 0; j < sizex; j++)
				{
					int val = *buf++ - 128;
					int index = val / 256;

					if (index >= 255)
						linei->Set(rcurve[index], gcurve[index], bcurve[index]);
					else
					{
						float sub = (val - index * 256) / 256.0f;
						linei->Set( rcurve[index] * (1-sub) + rcurve[index+1] * sub,
									gcurve[index] * (1-sub) + gcurve[index+1] * sub,
									bcurve[index] * (1-sub) + bcurve[index+1] * sub
								  );
									
					}

					linei++;
				}
			}
			else
			{
				for (j = 0; j < sizex; j++)
				{
					int val = buf[0] - 128;
					int index = val / 256;
					if (index >= 255)
						linei->r = rcurve[index];
					else
					{
						float sub = (val - index * 256) / 256.0f;
						linei->r = rcurve[index] * (1-sub) + rcurve[index+1] * sub;
					}

					val = buf[1] - 128;
					index = val / 256;
					if (index >= 255)
						linei->g = gcurve[index];
					else
					{
						float sub = (val - index * 256) / 256.0f;
						linei->g = gcurve[index] * (1-sub) + gcurve[index+1] * sub;
					}

					val = buf[2] - 128;
					index = val / 256;
					if (index >= 255)
						linei->b = rcurve[index];
					else
					{
						float sub = (val - index * 256) / 256.0f;
						linei->b = bcurve[index] * (1-sub) + bcurve[index+1] * sub;
					}

					linei++;
					buf += CHANNELS;
				} 
			}
		}
		_TIFFfree((void *) scanline);
	}  
	else if (BitsPerSample == 8)	// bytes
	{
		unsigned char *scanline = (unsigned char *) _TIFFmalloc(TIFFScanlineSize(tiffImage));
		_TIFFmemset(scanline, 0, TIFFScanlineSize(tiffImage));

		for (i = 0; i < sizey; i++)
		{
//			TIFF_Error_Flag = 0;
			TIFFReadScanline(tiffImage, scanline, i, 0);
//			if (TIFF_Error_Flag)
//				break;

			unsigned char *buf = scanline;
			fRGB *linei = line(i);

			if (CHANNELS==1)
			{
				for (j = 0; j < sizex; j++)
				{
					linei->Set(rcurve[*buf], gcurve[*buf], bcurve[*buf]);
					buf++;
					linei++;
				}
			}
			else
			{
				for (j = 0; j < sizex; j++)
				{
					linei->Set(rcurve[buf[0]], gcurve[buf[1]], bcurve[buf[2]]);
					linei++;
					buf += CHANNELS;
				} 
			}
		}
		_TIFFfree((void *) scanline);
	}  

	TIFFClose(tiffImage);
    
	sprintf(name, "%s", filename);

	return 1;
}


int floatimage::SaveTIF(const char *filename)
{
	TIFF *tiffImage = TIFFOpen(filename, "w");

	if (!tiffImage) 
		return 0;

	TIFFSetField(tiffImage, TIFFTAG_IMAGEWIDTH, sizex);
	TIFFSetField(tiffImage, TIFFTAG_IMAGELENGTH, sizey);
	TIFFSetField(tiffImage, TIFFTAG_ROWSPERSTRIP, sizey);
	TIFFSetField(tiffImage, TIFFTAG_SAMPLESPERPIXEL, 3);
	TIFFSetField(tiffImage, TIFFTAG_SAMPLEFORMAT, 3);
	TIFFSetField(tiffImage, TIFFTAG_COMPRESSION, 1);
	TIFFSetField(tiffImage, TIFFTAG_BITSPERSAMPLE, 32);
	TIFFSetField(tiffImage, TIFFTAG_PLANARCONFIG, 1);
	TIFFSetField(tiffImage, TIFFTAG_PHOTOMETRIC, 2);

		
	for (int row = 0; row < sizey; row++)
	{
		TIFFWriteRawStrip(tiffImage, 0, &line(row)[0], sizeof(float) * sizex * 3);
	}

	TIFFClose(tiffImage);

	return 1;
}


#endif

int floatimage::LoadBMP_Curve(const char *filename, float rcurve[256], float gcurve[256], float bcurve[256])
{
	FILE *infile = fopen(filename, "rb");

	if (infile == NULL)
		return 0;

	GFX_BITMAPINFOHEADER info;
	if (!ReadHeader(infile, &info))
	{
		fclose(infile);
		return 0;
	}

	Allocate(info.biWidth, info.biHeight);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	int stride = sizeof(RGB) * sizex;			// stride between lines in bytes
	stride = (stride + 3) & (~0x03);			// rounded up to nearest LONG boundary

	int filestride = (info.biBitCount / 8) * sizex;
	filestride = (filestride + 3) & (~0x03);

	if (info.biBitCount == 24)
	{
		RGB *buf = (RGB *) new char[stride];

		for (int y = 0; y < sizey; y++)
		{
			if (fread(buf, filestride, 1, infile) != 1)
			{
				Release();
				delete [] buf;
				fclose(infile);
				return 0;
			}

			RGB *cur = buf;
			fRGB *curline = line(sizey - 1 - y);
			for (int x = 0; x < sizex; x++)
			{
				curline[x].r = rcurve[cur->rgbRed];
				curline[x].g = gcurve[cur->rgbGreen];
				curline[x].b = bcurve[cur->rgbBlue];
				cur++;
			}
		}
		delete buf;
	}
	else if (info.biBitCount == 8)
	{
		unsigned char *buf = (unsigned char *) new char[filestride];

		for (int y = 0; y < sizey; y++)
		{
			if (fread(buf, filestride, 1, infile) != 1)
			{
				Release();
				delete [] buf;
				fclose(infile);
				return 0;
			}

			unsigned char *cur = buf;
			fRGB *curline = line(sizey - 1 - y);
			for (int x = 0; x < sizex; x++)
			{
				curline[x].r = rcurve[*cur];
				curline[x].g = gcurve[*cur];
				curline[x].b = bcurve[*cur];
				cur++;
			}
		}
		delete buf;
	}

	else if (info.biBitCount == 1)
	{
		unsigned char *buf = (unsigned char *) new unsigned char[sizex];

		for (int y = 0; y < sizey; y++)
		{
			if (fread(buf, sizex, 1, infile) != 1)
			{
				Release();
				delete [] buf;
				fclose(infile);
				return 0;
			}

			unsigned char *cur = buf;
			fRGB *curline = line(sizey - 1 - y);
			for (int x = 0; x < sizex; x++)
			{
				curline[x].r = *cur;
				curline[x].g = *cur;
				curline[x].b = *cur;
				cur++;
			}
		}
		delete buf;
	}
	
		//98366
		//98304
		//62 

	fclose(infile);

	sprintf(name, "%s", filename);

	return 1;
}

int floatimage::LoadBMP(const char *filename, float image_gamma)
{
	InitGamma(image_gamma);
	return LoadBMP_Curve(filename, gamma_table, gamma_table, gamma_table);
}

int LoadBMPChannel(buf2d<float> &fchannel, const char *filename)
{
	FILE *infile = fopen(filename, "rb");

	if (infile == NULL)
		return 0;

	GFX_BITMAPINFOHEADER info;
	if (!ReadHeader(infile, &info))
	{
		fclose(infile);
		return 0;
	}

	if (info.biBitCount != 8)
	{
		fclose(infile);
		return 0;
	}

	fchannel.Allocate(info.biWidth, info.biHeight);

	if (!fchannel.buf)
	{
		fclose(infile);
		return 0;
	}

	int stride = sizeof(RGB) * fchannel.sizex;	// stride between lines in bytes
	stride = (stride + 3) & (~0x03);			// rounded up to nearest LONG boundary

	int filestride = (info.biBitCount / 8) * fchannel.sizex;
	filestride = (filestride + 3) & (~0x03);
	int counter = 0;

	unsigned char *buf = (unsigned char *) new char[filestride];
	for (int y = 0; y < fchannel.sizey; y++)
	{
		counter++;
		if (fread(buf, filestride, 1, infile) != 1)
		{
			fchannel.Release();
			delete [] buf;
			fclose(infile);
			return 0;
		}

		unsigned char *cur = buf;
		float *curline = fchannel.line(fchannel.sizey - 1 - y);
		for (int x = 0; x < fchannel.sizex; x++)
		{
			curline[x] = *cur/255.0;
			cur++;
		}
	}
	delete buf;

	fclose(infile);
	return 1;
}

int LoadBMPLine(floatimage &fimg, int line, const char *filename)
{
	FILE *infile = fopen(filename, "rb");

	if (infile == NULL)
		return 0;

	GFX_BITMAPINFOHEADER info;
	if (!ReadHeader(infile, &info))
	{
		fclose(infile);
		return 0;
	}

	if(fimg.sizex != info.biWidth || fimg.sizey != 1)
		fimg.Allocate(info.biWidth, 1);

	if (!fimg.buf)
	{
		fclose(infile);
		return 0;
	}

	int stride = sizeof(RGB) * fimg.sizex;		// stride between lines in bytes
	stride = (stride + 3) & (~0x03);			// rounded up to nearest LONG boundary

	int filestride = (info.biBitCount / 8) * fimg.sizex;
	filestride = (filestride + 3) & (~0x03);

	if (info.biBitCount == 24)
	{
		RGB *buf = (RGB *) new char[stride];
		
		if(fseek(infile,(info.biHeight-line-1)*filestride,SEEK_CUR) != 0)
		{
			fimg.Release();
			delete [] buf;
			fclose(infile);
			return 0;
		}
		
		if (fread(buf, filestride, 1, infile) != 1)
		{
			fimg.Release();
			delete [] buf;
			fclose(infile);
			return 0;
		}

		RGB *cur = buf;
		fRGB *curline = fimg.line(0);
		for (int x = 0; x < fimg.sizex; x++)
		{
			curline[x].r = float(cur->rgbRed)/255.0;
			curline[x].g = float(cur->rgbGreen)/255.0;
			curline[x].b = float(cur->rgbBlue)/255.0;
			cur++;
		}
		
		delete buf;
	}

	fclose(infile);
	return 1;
}

/*
int floatimage::SaveBMP(const char *filename, float image_scale, float image_gamma)
{
	FILE *outfile = fopen(filename, "wb");

	GFX_BITMAPINFOHEADER info;
	info.biSize = sizeof(GFX_BITMAPINFOHEADER);
	info.biWidth = sizex;
	info.biHeight = sizey;
	info.biPlanes = 1;
	info.biBitCount = 24;
	info.biCompression = BI_RGB;
	info.biSizeImage = 0;
	info.biXPelsPerMeter = 2834;
	info.biYPelsPerMeter = 2834;
	info.biClrUsed = 0;
	info.biClrImportant = 0;

	if (!WriteHeader(&info, outfile))
	{
		fclose(outfile);
		return 0;
	}

	unsigned char lookup[65536];
	for (int x = 0; x < 65536; x++)
	{
		float tval = (float)(powf(short2float(x) * image_scale, 1.0 / image_gamma) * 256.f);
		lookup[x] = (tval > 255.0 ? 255 : (unsigned char) tval);
	}


	for (int y = sizey - 1; y >= 0; y--)
	{
		float *buf = &line(y)[0].r;
		for (x = 0; x < sizex; x++)
		{
			if ((fwrite(&lookup[float2short(buf[2])], 1, 1, outfile) != 1) ||
				(fwrite(&lookup[float2short(buf[1])], 1, 1, outfile) != 1) ||
				(fwrite(&lookup[float2short(buf[0])], 1, 1, outfile) != 1))
			{
				fclose(outfile);
				return 0;
			}

			buf += 3;
		}
	}
	
	fclose(outfile);
	return 1;
}
*/
int floatimage::SaveBMP(const char *filename, float image_scale, float image_gamma)
{
 FILE *outfile = fopen(filename, "wb");
 
 GFX_BITMAPINFOHEADER info;
 info.biSize = sizeof(GFX_BITMAPINFOHEADER);
 info.biWidth = sizex;
 info.biHeight = sizey;
 info.biPlanes = 1;
 info.biBitCount = 24;
 info.biCompression = BI_RGB;
 info.biSizeImage = 0;
 info.biXPelsPerMeter = 2834;
 info.biYPelsPerMeter = 2834;
 info.biClrUsed = 0;
 info.biClrImportant = 0;
 unsigned char blankspace[] = {0,0,0};
 
 int extra;
 
 
 if (!WriteHeader(&info, outfile))
 {
  fclose(outfile);
  return 0;
 }
 
 unsigned char lookup[65536];
 for (int x = 0; x < 65536; x++)
 {
  float tval = (float)(powf(short2float(x) * image_scale, 1.0 / image_gamma) * 256.f);
  lookup[x] = (tval > 255.0 ? 255 : (unsigned char) tval);
 }
 
 extra = sizex%4;
 for (int y = sizey - 1; y >= 0; y--)
 {
  float *buf = &line(y)[0].r;
  for (int x = 0; x < sizex; x++)
  {
   if ((fwrite(&lookup[float2short(buf[2])], 1, 1, outfile) != 1) ||
    (fwrite(&lookup[float2short(buf[1])], 1, 1, outfile) != 1) ||
    (fwrite(&lookup[float2short(buf[0])], 1, 1, outfile) != 1))
   {
    fclose(outfile);
    return 0;
   }
 
   buf += 3;
  }
  if (extra>0) fwrite(&blankspace,extra,1,outfile);
 }
 
 fclose(outfile);
 return 1;
}


int floatimage::SaveBMP_channel(const char *filename, int channel, float image_scale, float image_gamma)
{
	FILE *outfile = fopen(filename, "wb");

	GFX_BITMAPINFOHEADERWPALETTE info;

	info.info.biSize = sizeof(GFX_BITMAPINFOHEADER);
	info.info.biWidth = sizex;
	info.info.biHeight = sizey;
	info.info.biPlanes = 1;
	info.info.biBitCount = 8;
	info.info.biCompression = BI_RGB;
	info.info.biSizeImage = 0;
	info.info.biXPelsPerMeter = 2834;
	info.info.biYPelsPerMeter = 2834;
	info.info.biClrUsed = 256;
	info.info.biClrImportant = 256;

	for (int c = 0; c < 255; c++)
	{
		info.color[c].rgbRed = info.color[c].rgbGreen = info.color[c].rgbBlue = info.color[c].rgbReserved = c;
	}

	if (!WriteHeader(&info.info, outfile))
	{
		fclose(outfile);
		return 0;
	}

	unsigned char lookup[65536];
	for (int x = 0; x < 65536; x++)
	{
		float tval = (float)(powf(short2float(x) * image_scale, 1.0 / image_gamma) * 256.f);
		lookup[x] = (tval > 255.0 ? 255 : (unsigned char) tval);
	}

	for (int y = sizey - 1; y >= 0; y--)
	{
		float *buf = &line(y)[0].r;
		for (int x = 0; x < sizex; x++)
		{
			if (fwrite(&lookup[float2short(buf[channel])], 1, 1, outfile) != 1)
			{
				fclose(outfile);
				return 0;
			}

			buf += 3;
		}
	}
	
	fclose(outfile);
	return 1;
}

int floatimage::SavePPM(const char *filename, float image_scale, float image_gamma)
{
	FILE *outfile = fopen(filename, "wb");

	fputc('P', outfile);
	fputc('6', outfile);
	fputc(0x0a, outfile);

	fprintf(outfile, "%d %d", sizex, sizey);
	fputc(0x0a, outfile);

	fprintf(outfile, "255");
	fputc(0x0a, outfile);

	unsigned char lookup[65536];
	for (int x = 0; x < 65536; x++)
	{
		float tval = (float)(::pow((double) short2float(x) * image_scale,(double)  1.0 / image_gamma) * 256.f);
		lookup[x] = (tval > 255.0 ? 255 : (unsigned char) tval);
	}


	for (int y = 0; y < sizey; y++)
	{
		float *buf = &line(y)[0].r;
		for (int x = 0; x < sizex; x++)
		{
			if ((fwrite(&lookup[float2short(buf[0])], 1, 1, outfile) != 1) ||
				(fwrite(&lookup[float2short(buf[1])], 1, 1, outfile) != 1) ||
				(fwrite(&lookup[float2short(buf[2])], 1, 1, outfile) != 1))
			{
				fclose(outfile);
				return 0;
			}

			buf += 3;
		}
	}
	
	fclose(outfile);
	return 1;
}

int floatimage::SaveJPG(const char *filename, float image_scale, float image_gamma)
{
#ifndef FLOATIMAGE_USEJPEG
	return 0;
#else
	unsigned char lookup[65536];
	for (int x = 0; x < 65536; x++)
	{
		float tval = (float)(::pow((double) short2float(x) * image_scale,(double)  1.0 / image_gamma) * 256.f);
		lookup[x] = (tval > 255.0 ? 255 : (unsigned char) tval);
	}

	char *tmp = new char[sizex * sizey * 3];
	char *cur = tmp;

	if (!tmp)
		return 0;


	for (int y = 0; y < sizey; y++)
	{
		float *buf = &line(y)[0].r;
		for (x = 0; x < sizex; x++)
		{
			cur[0] = lookup[float2short(buf[0])];
			cur[1] = lookup[float2short(buf[1])];
			cur[2] = lookup[float2short(buf[2])];

			buf += 3; cur += 3;
		}
	}
	
	if(tmp == NULL)
		int per =1;
	int ok = JpegFile::RGBToJpegFile(filename, (unsigned char *)tmp, sizex, sizey, TRUE, 95);

	delete [] tmp;

	if (!ok) 
		return 0;

	return 1;
#endif
}

int floatimage::LoadPPM_Curve(const char *filename, float rcurve[256], float gcurve[256], float bcurve[256])
{
	FILE *infile = fopen(filename, "rb");

	if (infile == NULL)
		return 0;

	char temp[256];	

	if ((!fgets(temp, 256, infile)) || (temp[0] != 'P') || (temp[1] != '6'))
	{
		fclose(infile);
		return 0;
	}

	temp[0] = fgetc(infile);
	while (temp[0] == '#')
	{
		while (fgetc(infile) != '\n')
			;
//		fscanf(infile, "\n");
//		readline();
		temp[0] = fgetc(infile);
	}
	fgets(&temp[1], 255, infile);

	int width, height;
	sscanf(temp, "%d %d", &width, &height);

	int colors;
	fscanf(infile, "%d", &colors);
	fgetc(infile);

	if (colors != 255)
	{
		fclose(infile);
		return 0;
	}

	Allocate(width, height);
	if (!buf)
	{
		fclose(infile);
		return 0;
	}

	RGB *buf = new RGB[width];

	for (int y = 0; y < sizey; y++)
	{
		if (fread(buf, width * sizeof(RGB), 1, infile) != 1)
		{
			Release();
			delete [] buf;
			fclose(infile);
			return 0;
		}

		RGB *cur = buf;
		fRGB *curline = line(y);
		for (int x = 0; x < sizex; x++)
		{
			curline[x].r = rcurve[cur->rgbBlue];
			curline[x].g = gcurve[cur->rgbGreen];
			curline[x].b = bcurve[cur->rgbRed];
			cur++;
		}
	}

	delete buf;

	fclose(infile);

	sprintf(name, "%s", filename);

	return 1;
}


int floatimage::LoadPPM(const char *filename, float image_gamma)
{
	InitGamma(image_gamma);
	return LoadPPM_Curve(filename, gamma_table, gamma_table, gamma_table);
}


int floatimage::LoadJPG(const char *filename, float image_gamma)
{
#ifdef FLOATIMAGE_USEJPEG
	UINT height;
	UINT width;
	unsigned char *dataBuf;
	//read the file

	dataBuf = JpegFile::JpegFileToRGB(filename, &width, &height);
	if (dataBuf == NULL) 
		return 0;

	InitGamma(image_gamma);

	Allocate(width, height);

	unsigned char *cur = dataBuf;
	for (unsigned int y = 0; y < height; y++)
	{
		for (unsigned int x = 0; x < width; x++)
		{
			line(y)[x].Set(gamma_table[cur[0]], gamma_table[cur[1]], gamma_table[cur[2]]);
			cur += 3;
		}
	}

	delete [] dataBuf;
	return 1;
#else
	return 0;
#endif
}


int floatimage::LoadJPG_Curve(const char *filename, float rcurve[256], float gcurve[256], float bcurve[256])
{
#ifdef FLOATIMAGE_USEJPEG
	UINT height;
	UINT width;
	unsigned char *dataBuf;
	//read the file

	dataBuf = JpegFile::JpegFileToRGB(filename, &width, &height);
	if (dataBuf == NULL) 
		return 0;

	Allocate(width, height);

	unsigned char *cur = dataBuf;
	for (unsigned int y = 0; y < height; y++)
	{
		for (unsigned int x = 0; x < width; x++)
		{
			line(y)[x].Set(rcurve[cur[0]], gcurve[cur[1]], bcurve[cur[2]]);
			cur += 3;
		}
	}

	delete [] dataBuf;
	return 1;
#else
	return 0;
#endif
}

/*
void floatimage::SwapData(floatimage &in) 
{
	fRGB **tmp;
	int tmp_x,tmp_y;

	tmp = line; tmp_x = sizex; tmp_y = sizey;
	line = in.line; sizex = in.sizex; sizey = in.sizey;
	in.line = tmp; in.sizex = tmp_x; in.sizey = tmp_y;
}
*/
void floatimage::ResampleFrom(floatimage &in) 
{
	fRGB val;
	float *f = &line(0)[0].r;
	float inx,iny;
	int y,x;

	for (y = 0; y < sizey; y++) 
	{
		iny = ((y+0.5f)*in.sizey)/sizey-0.5f;
		for (x = 0; x < sizex; x++) 
		{
			inx = ((x+.5f)*in.sizex)/sizex-0.5f;

			val = in.GetPixelBilinear(inx,iny);

			*f++ = val.r;
			*f++ = val.g;
			*f++ = val.b;
		}
	}
}

void floatimage::gaussian_blur(int variance) {

	int i,count;
	float scalefac;
  
	scalefac = 1.f;
	count = 0;
   
	for (i=0;i<variance;i++) {
		blur1x();
		blur1y();
		count++;

		scalefac*=1.0f/16;

		if (count>5 || i==variance-1) {
			Scale(scalefac);
			count = 0; scalefac = 1.f;
		}

	}
}

void floatimage::fast_gaussian_blur(int variance, float *percent_complete) 
{
	int i,count;
	float variance_so_far,rough_target;
	float scalefac;
	int n;
  
	variance_so_far = 0.f;
	scalefac = 1.f;
	count = 0;

	n = 1;

	while (variance / m_boxvariance(n) > 15 ) n++;
	if (variance / m_boxvariance(n) < 5 ) n--;
	if (n==0) n=1;

//	printf("%f %d\n",m_boxvariance(n),n);

	float boxvariance = m_boxvariance(n);

	if (variance <= 3) rough_target = 0.f;
	else rough_target = (float)(variance - 4);

	if (rough_target<0.f) rough_target=0.f;

//	printf("target %f\n",rough_target);

	while (n > 1) 
	{
		int nextn = n / 2;
		float nextboxvariance = m_boxvariance(nextn);

		while (variance_so_far + boxvariance < rough_target - nextboxvariance) 
		{
			if (percent_complete)
				*percent_complete = variance_so_far / variance;

			this->box_blur(n);
//			printf("subblur %d\n",n);
			variance_so_far += boxvariance;
		}
		n = nextn;
		boxvariance = nextboxvariance;
	}
   
	for (i = (int) (variance_so_far); i < variance; i++) 
	{
		if (percent_complete)
			*percent_complete = i / (float) variance;

		blur1x();
		blur1y();
		count++;

//		printf("clean-up blur\n");

		scalefac*=1.0f/16;

		if (count > 5 || i == variance-1) 
		{
			Scale(scalefac);
			count = 0; scalefac = 1.f;
		}

	}
}


void floatimage::box_blur(int hw) {

	this->boxblurx(hw);
	this->boxblury(hw);
	this->Scale(1.f/((2*hw+1)*(2*hw+1)));

}

void floatimage::boxblurx(int hw) {
	int y,x;
	float *on,*left,*right;
	float *f = &line(0)[0].r;
	float *saveout,*out;

	saveout = (float*)malloc(sizex*3*sizeof(float));

	for (y = 0; y < sizey; y++) {

		left = f;
		on = f;
		right = f;
		out = saveout;

		out[0] = hw*left[0];
		out[1] = hw*left[1];
		out[2] = hw*left[2];

		for (x=0;x<=hw;x++) {
			out[0] += *right++;
			out[1] += *right++;
			out[2] += *right++;
		}
		out+=3;
		
		for (x=1;x<=hw;x++) {
			*out++ = *(out-3)+(*right++)-left[0];
			*out++ = *(out-3)+(*right++)-left[1];
			*out++ = *(out-3)+(*right++)-left[2];
		}
		for (x=hw+1;x<sizex-hw-1;x++) {
			*out++ = *(out-3)+(*right++)-(*left++);
			*out++ = *(out-3)+(*right++)-(*left++);
			*out++ = *(out-3)+(*right++)-(*left++);
		}
		for (x=sizex-hw-1;x<sizex;x++) {
			*out++ = *(out-3)+right[0]-(*left++);
			*out++ = *(out-3)+right[1]-(*left++);
			*out++ = *(out-3)+right[2]-(*left++);
		}

		memcpy(f,saveout,3*sizex*sizeof(float));

		f+=3*sizex;
	}

	free(saveout);
}

void floatimage::boxblury(int hw) {
	int y,x;
	float *on,*up,*down,*tmp;
	float *f = &line(0)[0].r;
	float *saveout,*out;

	saveout = (float*)malloc(sizey*3*sizeof(float));


	for (x = 0; x < sizex; x++) {

		up = f;
		on = f;
		down = f;
		out = saveout;

		out[0] = hw*up[0];
		out[1] = hw*up[1];
		out[2] = hw*up[2];

		for (y=0;y<=hw;y++) {
			out[0] += *down++;
			out[1] += *down++;
			out[2] += *down++;
			down+=3*(sizex-1);
		}
		out+=3;
		
		for (y=1;y<=hw;y++) {
			*out++ = *(out-3)+(*down++)-up[0];
			*out++ = *(out-3)+(*down++)-up[1];
			*out++ = *(out-3)+(*down++)-up[2];
			down+=3*(sizex-1);
		}
		for (y=hw+1;y<sizey-hw-1;y++) {
			*out++ = *(out-3)+(*down++)-(*up++);
			*out++ = *(out-3)+(*down++)-(*up++);
			*out++ = *(out-3)+(*down++)-(*up++);
			down+=3*(sizex-1); 	up+=3*(sizex-1);
		}
		for (y=sizey-hw-1;y<sizey;y++) {
			*out++ = *(out-3)+down[0]-(*up++);
			*out++ = *(out-3)+down[1]-(*up++);
			*out++ = *(out-3)+down[2]-(*up++);
			up+=3*(sizex-1);
		}
		tmp = f;
		for (y=0;y<sizey;y++) {
			*tmp++ = saveout[3*y+0];
			*tmp++ = saveout[3*y+1];
			*tmp++ = saveout[3*y+2];
			tmp+=3*(sizex-1);
		}

		f+=3;
	}

	free(saveout);

}

void floatimage::blur1x() {
  int  y;
  register float *on,*left,*right;
  register float *out;
  register int  x;
  float *f = &line(0)[0].r;
  float *saveout;

  saveout = (float*)malloc(sizex*3*sizeof(float));

  for (y = 0; y < sizey; y++) {
  
    out = saveout;
    on = f;
    left = f;
    right = f+3;
    *out++ = *on + *on + *left++ + *right++; on++;
    *out++ = *on + *on + *left++ + *right++; on++;
    *out++ = *on + *on + *left++ + *right++; on++;

    left-=3;
    for (x = 1; x < sizex-1; x++) {
		*out++ = *on + *on + *left++ + *right++; on++;
		*out++ = *on + *on + *left++ + *right++; on++;
		*out++ = *on + *on + *left++ + *right++; on++;
    }
	right-=3;
    *out++ = *on + *on + *left++ + *right++; on++;
    *out++ = *on + *on + *left++ + *right++; on++;
    *out++ = *on + *on + *left++ + *right++; on++;

	memcpy(f,saveout,3*sizex*sizeof(float));

	f+=3*sizex;

  }

  	free(saveout);

}

void floatimage::blur1y() {
  int  y;

  register float *on,*up,*down;
  register float *out;
  register int  x;
  float *f = &line(0)[0].r;
  float *saveout,*lastout,*tmp;

  saveout = (float*)malloc(sizex*3*sizeof(float));
  lastout = (float*)malloc(sizex*3*sizeof(float));

    out = lastout;
    on = f;
    up = f;
    down = f+3*sizex;
	for (x = 0; x < sizex; x++) {
		*out++ = *on + *on + *up++ + *down++; on++;
		*out++ = *on + *on + *up++ + *down++; on++;
		*out++ = *on + *on + *up++ + *down++; on++;
	}

    up-=3*sizex;
	tmp = saveout; saveout = lastout; lastout = tmp;

	for (y = 1; y < sizey-1; y++) {
		out = lastout;
		for (x = 0; x < sizex; x++) {
			*out++ = *on + *on + *up++ + *down++; on++;
			*out++ = *on + *on + *up++ + *down++; on++;
			*out++ = *on + *on + *up++ + *down++; on++;
		}
		memcpy(f,saveout,3*sizex*sizeof(float));
		tmp = saveout; saveout = lastout; lastout = tmp;
		f+=3*sizex;
	}

	down-=3*sizex;
	out = lastout;
	for (x = 0; x < sizex; x++) {
		*out++ = *on + *on + *up++ + *down++; on++;
		*out++ = *on + *on + *up++ + *down++; on++;
		*out++ = *on + *on + *up++ + *down++; on++;
	}

	memcpy(f,saveout,3*sizex*sizeof(float));
	memcpy(f+3*sizex,lastout,3*sizex*sizeof(float));

	free(lastout);
	free(saveout);



}

void floatimage::Scale(float cr,float cg,float cb) 
{
	int  y;
	register int  x;
	float *f = &line(0)[0].r;


	for (y = 0; y < sizey; y++) 
	{
		for (x = 0; x < sizex; x++) 
		{
			*f++ *= cr;
			*f++ *= cg;
			*f++ *= cb;
		}
	}

}

void floatimage::Scale(float s) 
{
	int  y;
	register int  x;
	float *f = &line(0)[0].r;

	for (y = 0; y < sizey; y++) 
	{
		for (x = 0; x < sizex; x++) 
		{
			*f++ *= s;
			*f++ *= s;
			*f++ *= s;
		}
	}
}

void floatimage::shoulder(float cut) 
{
	float scale,t,invscale,val;
	float *f = &line(0)[0].r;
	float *last = &line(sizey-1)[sizex-1].b;

	scale=2*(1-cut)/3.14159265f;
	invscale = 1/scale;

	while (f <= last)
	{
		val = *f;
		if (val>cut) 
		{
			t = (val-cut)*invscale;
			*f = scale*(float)atan(t)+cut;
		}
		f++;
	}
}

void floatimage::vignette(float r0, float r1) 
{
	int  y;
	register int  x;
	float radius;
	float reduce;
	float a,t;
	float *f = &line(0)[0].r;

	for (y = 0; y < sizey; y++) 
	{
		for (x = 0; x < sizex; x++) 
		{

			radius = (float)sqrt((double) (x-sizex/2)*(x-sizex/2) + (y-sizey/2)*(y-sizey/2));

			if (radius > r0) 
			{
				if (radius < r1) 
				{
					t = 1.0f - (radius-r0)/(r1-r0);
					a = t * 2 - 1;
					reduce = (float)((0.25 * M_PI + 0.5 * asin(a) + 0.5 * a * sqrt( 1 - a*a ))/(0.5 * M_PI));

				} else 
				{
					reduce =0.f;
				}
				*f++ *= reduce;
				*f++ *= reduce;
				*f++ *= reduce;
			}
			else 
			{
				f+=3;
			}
		}
	}
}


void floatimage::sepia(float cr,float cg,float cb) 
{
	int  x,y;
	float lum;

	for (y = 0; y < sizey; y++) 
	{
		for (x = 0; x < sizex; x++) 
		{

			lum = line(y)[x].Intensity();
			line(y)[x].Set(cr*lum,cg*lum,cb*lum);
		}
	}
}
/*
void floatimage::Rotate90CW(float *percent_complete)
{
	floatimage rotated;
	rotated.Allocate(sizey, sizex);

	if (rotated.line)
	{
		for (int y = 0; y < sizey; y++)
		{
			if (percent_complete)
				*percent_complete = (y / (float) sizey);
			for (int x = 0; x < sizex; x++)
				rotated.line(x)[sizey-1-y] = line(y)[x];
		}

		fRGB **temp = rotated.line;
		rotated.line = line;
		line = temp;

		sizey = rotated.sizey;
		sizex = rotated.sizex;
	}

	if (percent_complete)
		*percent_complete = 1.0f;
}
*/

void floatimage::ToSumTable()
{
	int x, y;

	for (x = 1; x < sizex; x++)
	{
		line(0)[x].r += line(0)[x-1].r;
		line(0)[x].g += line(0)[x-1].g;
		line(0)[x].b += line(0)[x-1].b;
	}

	for (y = 1; y < sizey; y++)
	{
		line(y)[0].r += line(y-1)[0].r;
		line(y)[0].g += line(y-1)[0].g;
		line(y)[0].b += line(y-1)[0].b;
	}

	for (y = 1; y < sizey; y++)
	{
		for (x = 1; x < sizex; x++)
		{
			line(y)[x].r += line(y-1)[x].r + line(y)[x-1].r - line(y-1)[x-1].r;
			line(y)[x].g += line(y-1)[x].g + line(y)[x-1].g - line(y-1)[x-1].g;
			line(y)[x].b += line(y-1)[x].b + line(y)[x-1].b - line(y-1)[x-1].b;
		}
	}
}


void floatimage::convolve_naive(floatimage &kernel, float *percent_complete)
{
	floatimage tmp;

	tmp.Allocate(sizex, sizey);

//	float	*rt,*gt,*bt;

	for (int y = 0; y < sizey; y++) 
	{
		if (percent_complete)
			*percent_complete = y / (float) sizey;
		for (int x = 0; x < sizex; x++) 
		{
			float	r,g,b;
			float	mr,mg,mb;

			r	=	0;
			g	=	0;
			b	=	0;
			mr	=	0;
			mb	=	0;
			mg	=	0;

			for (int ky = 0; ky < kernel.sizey; ky++) 
			{
				for (int kx = 0; kx < kernel.sizex; kx++)
				{
					int	cx = x + kx - kernel.sizex / 2;
					int	cy = y + ky - kernel.sizey / 2;

					if ((cx < 0) || (cy < 0) || (cx >= sizex) || (cy >= sizey)) 
					{
						mr += kernel.line(ky)[kx].r;
						mg += kernel.line(ky)[kx].g;
						mb += kernel.line(ky)[kx].b;
					}
					else 
					{
						r += kernel.line(ky)[kx].r * line(cy)[cx].r;
						g += kernel.line(ky)[kx].g * line(cy)[cx].g;
						b += kernel.line(ky)[kx].b * line(cy)[cx].b;
					}
				}
			}

			tmp.line(y)[x].r = r / (1.0f - mr);
			tmp.line(y)[x].g = g / (1.0f - mg);
			tmp.line(y)[x].b = b / (1.0f - mb);
		}
	}

	SwapData(tmp);

	if (percent_complete)
		*percent_complete = 1.0f;
}


void floatimage::convolve_naive_simple(floatimage &kernel)
{
	floatimage tmp;

	tmp.Allocate(sizex, sizey);

//	float	*rt,*gt,*bt;

	for (int y = 0; y < sizey; y++) 
	{
		for (int x = 0; x < sizex; x++) 
		{
			float	r,g,b;
//			float	mr,mg,mb;

			r	=	0;
			g	=	0;
			b	=	0;
//			mr	=	0;
//			mb	=	0;
//			mg	=	0;

			for (int ky = 0; ky < kernel.sizey; ky++) 
			{
				for (int kx = 0; kx < kernel.sizex; kx++)
				{
					int	cx = x + kx - kernel.sizex / 2;
					int	cy = y + ky - kernel.sizey / 2;

					if ((cx < 0) || (cy < 0) || (cx >= sizex) || (cy >= sizey)) 
					{
//						mr += kernel.line[ky][kx].r;
//						mg += kernel.line[ky][kx].g;
//						mb += kernel.line[ky][kx].b;
					}
					else 
					{
						r += kernel.line(ky)[kx].r * line(cy)[cx].r;
						g += kernel.line(ky)[kx].g * line(cy)[cx].g;
						b += kernel.line(ky)[kx].b * line(cy)[cx].b;
					}
				}
			}

			tmp.line(y)[x].r = r;
			tmp.line(y)[x].g = g;
			tmp.line(y)[x].b = b;
		}
	}

	SwapData(tmp);

}

void floatimage::UpsampleX(int newsizex)
{
/*	floatimage dest(newsizex, sizey);

	for (int x = 0; x < newsizex; x++)
	{
		float oldx = x * sizex / (newsizex);
	}
*/
}


void floatimage::DownsampleX(int newsizex)
{
}

void floatimage::BlitTo_NB(floatimage &dst, int dst_left, int dst_top, int src_left, int src_top, int src_width, int src_height)
{
	if (src_width < 0)
		src_width = sizex + src_width + 1;
	if (src_height < 0)
		src_height = sizey + src_height + 1;

	if (src_left < 0)
	{
		src_width += src_left;
		dst_left -= src_left;
		src_left = 0;
	}
	if (src_top < 0)
	{
		src_height += src_top;
		dst_top -= src_top;
		src_top = 0;
	}
	if (dst_left < 0)
	{
		src_left -= dst_left;
		src_width += dst_left;
		dst_left = 0;
	}
	if (dst_top < 0)
	{
		src_top -= dst_top;
		src_height += dst_top;
		dst_top = 0;
	}
	if (src_left + src_width > sizex)
	{
		src_width = sizex - src_left;
	}
	if (src_top + src_height > sizey)
	{
		src_height = sizey - src_top;
	}
	if (dst_left + src_width > dst.sizex)
	{
		src_width = dst.sizex - dst_left;
	}
	if (dst_top + src_height > dst.sizey)
	{
		src_height = dst.sizey - dst_top;
	}

	if (src_width && src_height)
		for (int y = 0; y < src_height; y++)
		{
			fRGB *curdst = &dst.line(dst_top + y)[dst_left];
			fRGB *cursrc = &line(src_top + y)[src_left];
			for (int x = 0; x < src_width; x++)
			{
				if (cursrc[x] != fRGB(0.0f, 0.0f, 0.0f))
					curdst[x] = cursrc[x];
			}
		}
}


void floatimage::AddTo(floatimage &dst, int dst_left, int dst_top, int src_left, int src_top, int src_width, int src_height)
{
	if (src_width < 0)
		src_width = sizex + src_width + 1;
	if (src_height < 0)
		src_height = sizey + src_height + 1;

	if (src_left < 0)
	{
		src_width += src_left;
		dst_left -= src_left;
		src_left = 0;
	}
	if (src_top < 0)
	{
		src_height += src_top;
		dst_top -= src_top;
		src_top = 0;
	}
	if (dst_left < 0)
	{
		src_left -= dst_left;
		src_width += dst_left;
		dst_left = 0;
	}
	if (dst_top < 0)
	{
		src_top -= dst_top;
		src_height += dst_top;
		dst_top = 0;
	}
	if (src_left + src_width > sizex)
	{
		src_width = sizex - src_left;
	}
	if (src_top + src_height > sizey)
	{
		src_height = sizey - src_top;
	}
	if (dst_left + src_width > dst.sizex)
	{
		src_width = dst.sizex - dst_left;
	}
	if (dst_top + src_height > dst.sizey)
	{
		src_height = dst.sizey - dst_top;
	}

	if (src_width && src_height)
		for (int y = 0; y < src_height; y++)
		{
			for (int x = 0; x < src_width; x++)
			{
				dst.line(dst_top + y)[dst_left + x] += line(src_top + y)[src_left + x];
			}
		}
}

#ifdef USE_IMAGEBMP

int floatimage::FromRGB(const imagergb &img, float image_gamma)
{
	gammatable gt(image_gamma);
	return FromRGB_Curve(img, gt.table, gt.table, gt.table);
}


int floatimage::FromRGB_Curve(const imagergb &img, float rcurve[256], float gcurve[256], float bcurve[256])
{
	Allocate(img.sizex, img.sizey);
	if (!buf)
	{
		return 0;
	}

	for (int y = 0; y < sizey; y++)
	{
		const bRGB *cur = img.line(y);
		fRGB *curline = line(y);
		for (int x = 0; x < sizex; x++)
		{
			curline[x].r = rcurve[cur->r];
			curline[x].g = gcurve[cur->g];
			curline[x].b = bcurve[cur->b];
			cur++;
		}
	}

	return 1;
}

#endif

fRGB floatimage::GetPixelBilinear(float sx, float sy)
{
	fRGB result;
	sx -= 0.49999f;
	sy -= 0.49999f;
	int bx = (int) sx;
	int by = (int) sy;
	int bxo = bx + 1;
	int byo = by + 1;
	float w00,w01,w10,w11;
	sx -= bx;
	sy -= by;

	if (bx<0) bx = bxo = 0;
	if (by<0) by = byo = 0;
	if (bx>=sizex-1) bx = bxo = sizex-1;
	if (by>=sizey-1) by = byo = sizey-1;

	w11 = sx*sy;

	w00 = 1 - sx - sy + w11;
	w10 = sx - w11;
	w01 = sy - w11;

	result.r = line(by)[bx].r * w00 +
			   line(by)[bxo].r * w10 +
			   line(byo)[bx].r * w01 +
			   line(byo)[bxo].r * w11;

	result.g = line(by)[bx].g * w00 +
			   line(by)[bxo].g * w10 +
			   line(byo)[bx].g * w01 +
			   line(byo)[bxo].g * w11;

	result.b = line(by)[bx].b * w00 +
			   line(by)[bxo].b * w10 +
			   line(byo)[bx].b * w01 +
			   line(byo)[bxo].b * w11;

	return result;
}

void floatimage::AddPixelBilinear(float sx, float sy, fRGB color)
{
	//fRGB result;
	sx -= 0.49999f;
	sy -= 0.49999f;
	int bx = (int) sx;
	int by = (int) sy;
	int bxo = bx + 1;
	int byo = by + 1;
	float w00,w01,w10,w11;
	sx -= bx;
	sy -= by;

	if (bx<0) bx = bxo = 0;
	if (by<0) by = byo = 0;
	if (bx>=sizex-1) bx = bxo = sizex-1;
	if (by>=sizey-1) by = byo = sizey-1;

	w11 = sx*sy;

	w00 = 1 - sx - sy + w11;
	w10 = sx - w11;
	w01 = sy - w11;

	line(by)[bx] += color*w00;
	line(by)[bxo] += color*w10;
	line(byo)[bx] += color*w01;
	line(byo)[bxo] += color*w11;
	
	/*
	result.r = line(by)[bx].r * w00 +
			   line(by)[bxo].r * w10 +
			   line(byo)[bx].r * w01 +
			   line(byo)[bxo].r * w11;

	result.g = line(by)[bx].g * w00 +
			   line(by)[bxo].g * w10 +
			   line(byo)[bx].g * w01 +
			   line(byo)[bxo].g * w11;

	result.b = line(by)[bx].b * w00 +
			   line(by)[bxo].b * w10 +
			   line(byo)[bx].b * w01 +
			   line(byo)[bxo].b * w11;

	return result;
	*/
}

floatimage floatimage::HalfSizeDeleteOld()
{
	floatimage temp;
	temp.Allocate(sizex / 2, sizey / 2);

	for (int y = 0; y < temp.sizey; y++)
	{
		int cy = 2 * y;
		for (int x = 0; x < temp.sizex; x++)
		{
			int cx = 2 * x;
			temp.line(y)[x] = (line(cy)[cx] + line(cy+1)[cx] + line(cy)[cx+1] + line(cy+1)[cx+1])*0.25f;
		}
	}
	delete this;
	return temp;
}

void floatimage::HalfSize()
{
	floatimage temp;
	temp.Allocate(sizex / 2, sizey / 2);

	for (int y = 0; y < temp.sizey; y++)
	{
		int cy = 2 * y;
		for (int x = 0; x < temp.sizex; x++)
		{
			int cx = 2 * x;
			temp.line(y)[x] = (line(cy)[cx] + line(cy+1)[cx] + line(cy)[cx+1] + line(cy+1)[cx+1])*0.25f;
		}
	}

	SwapData(temp);
}

void floatimage::DoubleSize()
{
	floatimage temp;
	temp.Allocate(sizex * 2, sizey * 2);

	for (int y = 0; y < temp.sizey; y++)
	{
		float cy = (y + 0.5f) / 2.f;
		for (int x = 0; x < temp.sizex; x++)
		{
			float cx = (x + 0.5f) / 2.f;
			temp.line(y)[x] = GetPixelBilinear(cx,cy);
		}
	}

	SwapData(temp);
}

void floatimage::Clamp(float r0, float g0, float b0, float r1, float g1, float b1) {
	int y,x;
	for (y=0;y<sizey;y++) {
		for (x=sizex;x>=0;x--) {
			if (line(y)[x].r > r1) line(y)[x].r = r1;
			if (line(y)[x].g > g1) line(y)[x].g = g1;
			if (line(y)[x].b > b1) line(y)[x].b = b1;
			if (line(y)[x].r < r0) line(y)[x].r = r0;
			if (line(y)[x].g < g0) line(y)[x].g = g0;
			if (line(y)[x].b < b0) line(y)[x].b = b0;

		}
	}
}

/*void floatimage::Clear(float r0, float g0, float b0) {
	int y,x;
	for (y=0;y<sizey;y++) {
		for (x=0;x<sizex;x++) {
			line(y)[x].r = r0;
			line(y)[x].g = g0;
			line(y)[x].b = b0;
		}
	}
}*/

void floatimage::Quantize(int levels, float minval, float maxval)
{
	float diff = maxval - minval;
	float *f = &line(0)[0].r;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		if (*f < minval)
			*f = minval;
		else if (*f > maxval)
			*f = maxval;
		else
			*f = ((int)((*f - minval) * levels / diff + 0.5)) * diff / levels + minval;
		f++;
	}			
}

int floatimage::CopyChannel(floatimage &other, int src_channel, int dst_channel)
{
	if ((sizex != other.sizex) || (sizey != other.sizey))
		return 0;
	
	float *f = line(0)[0].channel(dst_channel);
	float *o = other.line(0)[0].channel(src_channel);
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		*f++ = *o++;
		f++; o++;
		f++; o++;
	}
	
	return 1;
}



#ifdef FLOATIMAGE_USEIMAGE
void floatimage::addWithGamma(image &i, float gamma = 1.75f)
{
	InitGamma(gamma);
	
	float *f = &line(0)[0].r;
	unsigned char *b = &i.line(0)[0].rgbBlue;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		f[2] += gamma_table[*b++];
		f[1] += gamma_table[*b++];
		f[0] += gamma_table[*b++];
		f += 3;
	}
}

void floatimage::minWithGamma(image &i, float gamma = 1.75f)
{
	InitGamma(gamma);
	
	float *f = &line(0)[0].r;
	unsigned char *b = &i.line(0)[0].rgbBlue;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		f[2] = min(f[2], gamma_table[*b]);  b++;
		f[1] = min(f[1], gamma_table[*b]);  b++;
		f[0] = min(f[0], gamma_table[*b]);  b++;
		f += 3;
	}
}

void floatimage::ToImage(image &i, float gamma = 1.75f)
{
	InitGamma(gamma);
	
	i.Allocate(sizex, sizey);
	
	float *f = &line(0)[0].r;
	unsigned char *b = &i.line(0)[0].rgbBlue;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		b[0] = f[2] * 255.0f; // blue
		b[1] = f[1] * 255.0f; // green
		b[2] = f[0] * 255.0f; // red
		
		b += 3;
		f += 3;
	}
}

#endif

/*	operator *= (float s)
{
float *f = &line(0)[0].r;
float *last = &line(sizey-1)[sizex-1].b;
while (f <= last)
{
*f++ *= s;
*f++ *= s;
*f++ *= s;
}
}

		operator -= (float s)
		{
		float *f = &line(0)[0].r;
		float *last = &line(sizey-1)[sizex-1].b;
		while (f <= last)
		{
		*f++ -= s;
		*f++ -= s;
		*f++ -= s;
		}
		}*/

void floatimage::Blend(floatimage &ia, floatimage &ib, floatimage &ic, floatimage &id)
{
	Allocate(ia.sizex, ia.sizey);
	
	float *d = &line(0)[0].r;
	float *sa = &ia.line(0)[0].r;
	float *sb = &ib.line(0)[0].r;
	float *sc = &ic.line(0)[0].r;
	float *sd = &id.line(0)[0].r;
	float *dlast = &line(sizey-1)[sizex-1].b;
	while (d <= dlast)
	{
		int count = 0;
		d[0] = d[1] = d[2] = 0.f;
		
		if (sa[0] > 0.0f)
		{
			d[0] += sa[0]; d[1] += sa[1]; d[2] += sa[2];
			count++;
		}
		
		if (sb[0] > 0.0f)
		{
			d[0] += sb[0]; d[1] += sb[1]; d[2] += sb[2];
			count++;
		}
		
		if (sc[0] > 0.0f)
		{
			d[0] += sc[0]; d[1] += sc[1]; d[2] += sc[2];
			count++;
		}
		
		if (sd[0] > 0.0f)
		{
			d[0] += sd[0]; d[1] += sd[1]; d[2] += sd[2];
			count++;
		}
		
		if (count)
		{
			d[0] /= count;
			d[1] /= count;
			d[2] /= count;
		}
		else
		{
			d[0] = -FLT_MAX;
			d[1] = -FLT_MAX;
			d[2] = -FLT_MAX;
		}
		
		d += 3; sa += 3; sb += 3; sc += 3; sd += 3;
	}			
}


void floatimage::pow(double p)
{
	float *f = &line(0)[0].r;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		f[0] = (float)::pow((double) f[0],(double) p);
		f[1] = (float)::pow((double) f[1],(double) p);
		f[2] = (float)::pow((double) f[2],(double) p);
		f += 3;
	}
}


void floatimage::MultiplyInverse(floatimage &a, floatimage &b)
{
	Allocate(a.sizex, a.sizey);
	
	float *d = &line(0)[0].r;
	float *sa = &a.line(0)[0].r;
	float *sb = &b.line(0)[0].r;
	float *dlast = &line(sizey-1)[sizex-1].b;
	while (d <= dlast)
	{
		*d++ = *sa++ * (1.0f - *sb++);		// WOOOOOOOOOOOOOOOOOOOOOOOO!!!!!!!!!!!!!
		*d++ = *sa++ * (1.0f - *sb++);
		*d++ = *sa++ * (1.0f - *sb++);
	}			
}


void floatimage::absolutevalue()
{
	float *f = &line(0)[0].r;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		*f = (float) fabs(*f); f++;
		*f = (float) fabs(*f); f++;
		*f = (float) fabs(*f); f++;
	}
}

void floatimage::ClampRange(float min, float max)
{
	float *f = &line(0)[0].r;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		if (*f < min)
			*f = min;
		else if (*f > max)
			*f = max;
		f++;
		if (*f < min)
			*f = min;
		else if (*f > max)
			*f = max;
		f++;
		if (*f < min)
			*f = min;
		else if (*f > max)
			*f = max;
		f++;
	}
}

void floatimage::minimum(floatimage &other)
{
	float *f = &line(0)[0].r;
	float *o = &other.line(0)[0].r;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		f[0] = min(f[0], o[0]);
		f[1] = min(f[1], o[1]);
		f[2] = min(f[2], o[2]);
		f += 3; o += 3;
	}
}

void floatimage::maximum(floatimage &other)
{
	float *f = &line(0)[0].r;
	float *o = &other.line(0)[0].r;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		f[0] = max(f[0], o[0]);
		f[1] = max(f[1], o[1]);
		f[2] = max(f[2], o[2]);
		f += 3; o += 3;
	}
}

void floatimage::find_maximum(int x0, int y0, int x1, int y1, int &maxx, int &maxy)
{
	float damax = line(y0)[x0].Intensity();
	maxx = x0; maxy = y0;
	
	for (int y = y0; y <= y1; y++)
	{
		fRGB *curline = line(y);
		for (int x = x0; x <= x1; x++)
		{
			float intensity = curline[x].Intensity();
			if (intensity > damax)
			{
				damax = intensity;
				maxx = x;
				maxy = y;
			}
		}
	}
}


void floatimage::find_minimum(int x0, int y0, int x1, int y1, int &minx, int &miny)
{
	float damin = line(y0)[x0].Intensity();
	minx = x0; miny = y0;
	
	for (int y = y0; y <= y1; y++)
	{
		fRGB *curline = line(y);
		for (int x = x0; x <= x1; x++)
		{
			float intensity = curline[x].Intensity();
			if (intensity < damin)
			{
				damin = intensity;
				minx = x;
				miny = y;
			}
		}
	}
}

//		float fit_quadratic_max(float xn1, float x0, float x1)
//		{
//			float a = (xn1 + x1) / 2 - x0;
//			float b = (x1 - xn1) / 2;
//			float c = x0;
//			return -b / (2 * a);
//
//			return (xn1 - x1) / (2 * (xn1 + x1) - 4 * x0);
//		}

void floatimage::find_maximum_subpixel(float &maxx, float &maxy)
{
	//			while change
	//				find_maximum(maxx-1, maxy-1, maxx+1, maxy+1, maxx, maxy);
	
	int x = (int) maxx;
	int y = (int) maxy;
	
	if ((x <= 0) || (x >= sizex - 1) || (y <= 0) || (y >= sizey-1))
		return;
	
		/*			float sx = maxx - x;
		float sy = maxy - y;
		
		  int x0, x1, y0, y1;
		  
			if (sx < 0.5)
			x0 = x - 1;
			else
			x0 = x;
			
			  if (sy < 0.5)
			  y0 = y - 1;
			  else
			  y0 = y;
			  x1 = x0 + 1; y1 = y0 + 1;
	*/
	
#define find_quadratic_minmax(xn1, x0, x1)    (((xn1) - (x1)) / (2 * ((xn1) + (x1)) - 4 * (x0)))
	
	float vn1 = line(y)[x-1].Intensity();
	float v0  = line(y)[x  ].Intensity();
	float v1  = line(y)[x+1].Intensity();
	
	float nsx = find_quadratic_minmax(vn1, v0, v1) + 0.5f;
	
	vn1 = line(y-1)[x].Intensity();
	v0  = line(y  )[x].Intensity();
	v1  = line(y+1)[x].Intensity();
	
	float nsy = find_quadratic_minmax(vn1, v0, v1) + 0.5f;
	
	maxx = x + nsx;
	maxy = y + nsy;
	
#undef find_quadratic_minmax
}

fRGB floatimage::maximum(int x0, int y0, int x1, int y1)
{
	fRGB damax = line(y0)[x0];
	
	for (int y = y0; y <= y1; y++)
	{
		for (int x = x0; x <= x1; x++)
		{
			if (line(y)[x].r > damax.r)
				damax.r = line(y)[x].r;
			if (line(y)[x].g > damax.g)
				damax.g = line(y)[x].g;
			if (line(y)[x].b > damax.b)
				damax.b = line(y)[x].b;
		}
	}
	
	return damax;
}


fRGB floatimage::minimum(int x0, int y0, int x1, int y1)
{
	fRGB damin = line(y0)[x0];
	
	for (int y = y0; y <= y1; y++)
	{
		for (int x = x0; x <= x1; x++)
		{
			if (line(y)[x].r < damin.r)
				damin.r = line(y)[x].r;
			if (line(y)[x].g < damin.g)
				damin.g = line(y)[x].g;
			if (line(y)[x].b < damin.b)
				damin.b = line(y)[x].b;
		}
	}
	
	return damin;
}

void floatimage::ScaleToRange(float min, float max)
{
	//			fRGB mincolors = minimum(0, 0, sizex-1, sizey-1);
	fRGB maxcolors = maximum(0, 0, sizex-1, sizey-1);
	//			float tmin = min(min(mincolors.r, mincolors.g), mincolors.b);
	float tmax = max(max(maxcolors.r, maxcolors.g), maxcolors.b);
	
	*this *= (max / tmax);
}

fRGB floatimage::average(int x0, int y0, int x1, int y1)
{
	int count = 0;
	fRGB dasum(0.0f, 0.0f, 0.0f);
	
	for (int y = y0; y <= y1; y++)
	{
		for (int x = x0; x <= x1; x++)
		{
			dasum.r += line(y)[x].r;
			dasum.g += line(y)[x].g;
			dasum.b += line(y)[x].b;
			count++;
		}
	}
	
	dasum.r /= count;
	dasum.g /= count;
	dasum.b /= count;
	
	return dasum;
}


fRGB floatimage::SumTable_sum(int x0, int y0, int x1, int y1)
{
	fRGB dasum;
	
	dasum = line(y1)[x1];
	
	if (x0 > 0)
		dasum -= line(y1)[x0-1];
	
	if (y0 > 0)
	{
		dasum -= line(y0-1)[x1];
		
		if (x0 > 0)		// both - must subtract overlap
		{
			dasum += line(y0-1)[x0-1];
		}
	}
	
	return dasum;
}

fRGB floatimage::SumTable_average(int x0, int y0, int x1, int y1)
{
	fRGB dasum;
	
	dasum = line(y1)[x1];
	
	if (x0 > 0)
		dasum -= line(y1)[x0-1];
	
	if (y0 > 0)
	{
		dasum -= line(y0-1)[x1];
		
		if (x0 > 0)		// both - must subtract overlap
		{
			dasum += line(y0-1)[x0-1];
		}
	}
	
	return dasum / (float) ((x1 - x0 + 1) * (y1 - y0 + 1));
}

fRGB floatimage::GetPixelBilinearSlowClip(float sx, float sy)
{
	fRGB result;
	int bx = (int) sx;
	int by = (int) sy;
	int bxo = bx + 1;
	int byo = by + 1;
	sx -= bx;
	sy -= by;
	
	if (bx >= sizex)
		bx = sizex - 1;
	if (by >= sizey)
		by = sizey - 1;
	if (bxo >= sizex)
		bxo = sizex - 1;
	if (byo >= sizey)
		byo = sizey - 1;
	if (bx < 0)
		bx = 0;
	if (by < 0)
		by = 0;
	if (bxo < 0)
		bxo = 0;
	if (byo < 0)
		byo = 0;
	
	float osx = 1.0f - sx;
	float osy = 1.0f - sy;
	
	result.r = line(by)[bx].r * osx * osy +
					   line(by)[bxo].r * sx * osy +
					   line(byo)[bx].r * osx * sy +
					   line(byo)[bxo].r * sx * sy;
	
	result.g = line(by)[bx].g * osx * osy +
					   line(by)[bxo].g * sx * osy +
					   line(byo)[bx].g * osx * sy +
					   line(byo)[bxo].g * sx * sy;
	
	result.b = line(by)[bx].b * osx * osy +
					   line(by)[bxo].b * sx * osy +
					   line(byo)[bx].b * osx * sy +
					   line(byo)[bxo].b * sx * sy;
	
	return result;
}


fRGB floatimage::GetPixelBilinearLogSlow(float sx, float sy)		// cubic, not energy-preserving, non-negative, interpolating
{
	fRGB result;
	
	if (sx < 0.0f)
		sx = 0.0f;
	else if (sx >= sizex - 1.0f)
		sx = sizex - 1.0f;
	if (sy < 0.0f)
		sy = 0.0f;
	else if (sy >= sizey - 1.0f)
		sy = sizey - 1.0f;
	
	int bx = (int) sx;
	int by = (int) sy;
	int bxo = bx + 1;  if (bxo > sizex - 1)  bxo = sizex - 1;
	int byo = by + 1;  if (byo > sizey - 1)  byo = sizey - 1;
	
	sx -= bx;
	sy -= by;
	
	float osx = 1.0f - sx;
	float osy = 1.0f - sy;
	
	result.r = (float)exp(log(line(by)[bx].r + 0.0001) * osx * osy +
		log(line(by)[bxo].r + 0.0001) * sx * osy +
		log(line(byo)[bx].r + 0.0001) * osx * sy +
		log(line(byo)[bxo].r + 0.0001) * sx * sy) - 0.0001f;
	
	result.g = (float)exp(log(line(by)[bx].g + 0.0001) * osx * osy +
		log(line(by)[bxo].g + 0.0001) * sx * osy +
		log(line(byo)[bx].g + 0.0001) * osx * sy +
		log(line(byo)[bxo].g + 0.0001) * sx * sy) - 0.0001f;
	
	result.b = (float)exp(log(line(by)[bx].b + 0.0001) * osx * osy +
		log(line(by)[bxo].b + 0.0001) * sx * osy +
		log(line(byo)[bx].b + 0.0001) * osx * sy +
		log(line(byo)[bxo].b + 0.0001) * sx * sy) - 0.0001f;
	
	return result;
}

fRGB floatimage::GetPixelBSplineSlow(float sx, float sy)	// cubic, energy-preserving, non-negative, non-interpolating, slightly 'blurry'
{
	float A = -0.75f;
	fRGB val(0, 0, 0);
	float totalweight = 0;
	
	if (sx < 0.0f)
		sx = 0.0f;
	else if (sx >= sizex - 1.0f)
		sx = sizex - 1.0f;
	if (sy < 0.0f)
		sy = 0.0f;
	else if (sy >= sizey - 1.0f)
		sy = sizey - 1.0f;
	
	int bx = (int) sx;
	int by = (int) sy;
	
	//			bx %= sizex;
	//			by %= sizey;
	
	for (int y = by - 1; y < by + 3; y++)
	{
		if ((y >= 0) && (y < sizey))
		{
			float disty = (float)fabs(sy - y);
			float yweight;
			
			if (disty <= 1)			// b-spline
				yweight = (4.f + (-6.f + 3.f * disty) * disty * disty) / 6.f;
			else
				yweight = (2.f - disty) * (2.f - disty) * (2.f - disty) / 6.f;
			
			for (int x = bx - 1; x < bx + 3; x++)
			{
				if ((x >= 0) && (x < sizex))
				{
					float distx = (float)fabs(sx - x);
					float weight;
					
					if (distx <= 1)			// b-spline
						weight = ((4.f + (-6.f + 3.f * distx) * distx * distx) / 6.f) * yweight;
					else
						weight = ((2.f - distx) * (2.f - distx) * (2.f - distx) / 6.f) * yweight;
					
					val += line(y)[x] * weight;
					totalweight += weight;
				}
			}
		}
	}
	val *= 1.0f / totalweight;
	return val;
}


fRGB floatimage::GetPixelBicubicLogSlow(float sx, float sy)
{
	float A = -0.75f;
	fRGB val(0, 0, 0);
	float totalweight = 0;
	
	if (sx < 0.0f)
		sx = 0.0f;
	else if (sx >= sizex - 1.0f)
		sx = sizex - 1.0f;
	if (sy < 0.0f)
		sy = 0.0f;
	else if (sy >= sizey - 1.0f)
		sy = sizey - 1.0f;
	
	int bx = (int) sx;
	int by = (int) sy;
	
	//			bx %= sizex;
	//			by %= sizey;
	
	for (int y = by - 1; y < by + 3; y++)
	{
		if ((y >= 0) && (y < sizey))
		{
			float disty = (float)fabs(sy - y);
			float yweight;
			
			if (disty <= 1)			// photoshop bicubic
				yweight = (( A + 2.f ) * disty - ( A + 3.f )) * disty * disty + 1.f;
			else
				yweight = (( A * disty - 5.f * A ) * disty + 8.f * A ) * disty - 4.f * A;
			
				/*
				if (disty <= 1)			// 'true' bicubic ?
				yweight  = ( ( disty - 9.0/5.0 ) * disty -   1.0/5.0     ) * disty + 1.0;
				else
				yweight  = ( ( -1.0/3.0 * (disty-1) + 4.0/5.0     ) * (disty-1) -   7.0/15.0 ) * (disty-1);
			*/
			for (int x = bx - 1; x < bx + 3; x++)
			{
				if ((x >= 0) && (x < sizex))
				{
					float distx = (float)fabs(sx - x);
					float weight;
					
					if (distx <= 1)		// photoshop bicubic
						weight = ((( A + 2.f ) * distx - (A + 3.f)) * distx * distx +1.f) * yweight;
					else
						weight = ((( A * distx - 5.f * A ) * distx + 8.f * A ) * distx - 4.f * A) * yweight;
					
						/*							if (distx <= 1)		// 'true' bicubic ?
						weight  = (( ( distx - 9.0/5.0 ) * distx -   1.0/5.0     ) * distx + 1.0) * yweight;
						else
						weight  = (( ( -1.0/3.0 * (distx-1) + 4.0/5.0     ) * (distx-1) -   7.0/15.0 ) * (distx-1)) * yweight;
					*/
					val.r += (float)log(line(y)[x].r + 0.0001) * weight;
					val.g += (float)log(line(y)[x].g + 0.0001) * weight;
					val.b += (float)log(line(y)[x].b + 0.0001) * weight;
					//							val += line(y)[x] * weight;
					totalweight += weight;
				}
			}
		}
	}
	val *= 1.0f / totalweight;
	val.r = (float)(exp(val.r) - 0.0001f);
	val.g = (float)(exp(val.g) - 0.0001f);
	val.b = (float)(exp(val.b) - 0.0001f);
	return val;
}


void floatimage::ResampleX(int newsizex)
{
	if (newsizex > sizex)
		UpsampleX(newsizex);
	else
		DownsampleX(newsizex);
}

void floatimage::GetSubImage(int xstart, int ystart, int xwidth, int yheight, floatimage &dest)
{
	dest.Allocate(xwidth, yheight);
	for (int y = 0; y < yheight; y++)
		for (int x = 0; x < xwidth; x++)
		{
			dest.line(y)[x] = line(ystart + y)[xstart + x];
		}
}

void floatimage::AddBlitModulate(int xstart, int ystart, floatimage &source, float r, float g, float b)
{
	for (int y = 0; y < source.sizey; y++)
	{
		int yval = ystart + y;
		if ((yval < 0) || (yval >= sizey))
			continue;
		for (int x = 0; x < source.sizex; x++)
		{
			int xval = xstart + x;
			if ((xval < 0) || (xval >= sizex))
				continue;
			
			line(yval)[xval].r += source.line(y)[x].r * r;
			line(yval)[xval].g += source.line(y)[x].g * g;
			line(yval)[xval].b += source.line(y)[x].b * b;
		}
	}
}

void floatimage::AddBlitScaleModulate(int xstart, int ystart, float scale, floatimage &source, float r, float g, float b)
{
	int width = (int) (source.sizex * scale);
	int height = (int) (source.sizey * scale);
	
	xstart -= width / 2;
	ystart -= height / 2;
	
	floatimage blah;
	
	if (scale < 1.0)
	{
		blah = source;
		source.SwapData(blah);
		source.fast_gaussian_blur((int) ((0.7 / scale) * (0.7 / scale)));
		//				source.box_blur(0.5 / scale);
	}
	
	for (int y = 0; y < height; y++)
	{
		int yval = ystart + y;
		if ((yval < 0) || (yval >= sizey))
			continue;
		for (int x = 0; x < width; x++)
		{
			int xval = xstart + x;
			if ((xval < 0) || (xval >= sizex))
				continue;
			
			fRGB temp = source.GetPixelBilinearSlow(x / scale, y / scale);
			line(yval)[xval].r += temp.r * r;
			line(yval)[xval].g += temp.g * g;
			line(yval)[xval].b += temp.b * b;
		}
	}
	
	if (scale < 1.0)
	{
		source.SwapData(blah);
		blah.Release();
	}
}

void floatimage::AddBlitScaleFastSubPixel(float xstart, float ystart, float scale, floatimage &source)
{
	float width = source.sizex * scale;
	float height = source.sizey * scale;
	
	xstart -= width / 2;
	ystart -= height / 2;
	
	for (int y = 0; y < height; y++)
	{
		int yval = (int) (ystart + y);
		if ((yval < 0) || (yval >= sizey))
			continue;
		for (int x = 0; x < width; x++)
		{
			int xval = (int) (xstart + x);
			if ((xval < 0) || (xval >= sizex))
				continue;
			
			fRGB temp = source.line((int) (y / scale))[(int) (x / scale)];
			line(yval)[xval].r += temp.r;
			line(yval)[xval].g += temp.g;
			line(yval)[xval].b += temp.b;
		}
	}
}

void floatimage::AddBlitScale(int xstart, int ystart, float scale, floatimage &source)
{
	int width = (int) (source.sizex * scale);
	int height = (int) (source.sizey * scale);
	
	xstart -= width / 2;
	ystart -= height / 2;
	
	for (int y = 0; y < height; y++)
	{
		int yval = ystart + y;
		if ((yval < 0) || (yval >= sizey))
			continue;
		for (int x = 0; x < width; x++)
		{
			int xval = xstart + x;
			if ((xval < 0) || (xval >= sizex))
				continue;
			
			fRGB temp = source.line((int) (y / scale))[(int) (x / scale)];
			line(yval)[xval].r += temp.r;
			line(yval)[xval].g += temp.g;
			line(yval)[xval].b += temp.b;
		}
	}
}

void floatimage::AddBlitScaleFast(int xstart, int ystart, float scale, floatimage &source)
{
	int width = (int) (source.sizex * scale);
	int height = (int) (source.sizey * scale);
	
	xstart -= width / 2;
	ystart -= height / 2;
	
	floatimage blah;
	
	if (scale < 1.0)
	{
		blah = source;
		source.SwapData(blah);
		//				source.fast_gaussian_blur((0.7 / scale) * (0.7 / scale));
		source.box_blur((int) (3.0 / scale));
	}
	
	for (int y = 0; y < height; y++)
	{
		int yval = ystart + y;
		if ((yval < 0) || (yval >= sizey))
			continue;
		for (int x = 0; x < width; x++)
		{
			int xval = xstart + x;
			if ((xval < 0) || (xval >= sizex))
				continue;
			
			fRGB temp = source.GetPixelBilinearSlow(x / scale, y / scale);
			line(yval)[xval].r += temp.r;
			line(yval)[xval].g += temp.g;
			line(yval)[xval].b += temp.b;
		}
	}
	
	if (scale < 1.0)
	{
		source.SwapData(blah);
		blah.Release();
	}
}

void floatimage::ShiftX(float percent, floatimage &dest)
{
	while (percent < 0.0)
		percent += 1.0;
	float shift = percent * sizex;
	dest.Allocate(sizex, sizey);
	for (int y = 0; y < sizey; y++)
		for (int x = 0; x < sizex; x++)
		{
			float sx = x + shift;
			
			//					while (sx < 0.0)
			//						sx += sizex;
			//					while (sx >= sizex)
			//						sx -= sizex;
			
			if (sx < 0.0)
				sx += ((int) (sx / sizex) + 1) * sizex;
			
			sx -= ((int) (sx / sizex)) * sizex;
			
			dest.line(y)[x] = GetPixelBilinearSlow(sx, y);
		}
}

int floatimage::CountColor(float r, float g, float b)
{
	int count = 0;
	float *f = &line(0)[0].r;
	float *last = &line(sizey-1)[sizex-1].b;
	while (f <= last)
	{
		if ((f[0] == r) && (f[1] == g) && (f[2] == b))
		{
			count++;
		}
		f += 3;
	}
	
	return count;
}

void floatimage::AddPixelBilinearSlowClip(float sx, float sy,fRGB c)
{
	int bx = (int) sx;
	int by = (int) sy;
	int bxo = bx + 1;
	int byo = by + 1;
	sx -= bx;
	sy -= by;
	
	if (bx >= sizex)
		bx = sizex - 1;
	if (by >= sizey)
		by = sizey - 1;
	if (bxo >= sizex)
		bxo = sizex - 1;
	if (byo >= sizey)
		byo = sizey - 1;
	if (bx < 0)
		bx = 0;
	if (by < 0)
		by = 0;
	if (bxo < 0)
		bxo = 0;
	if (byo < 0)
		byo = 0;
	
	float osx = 1.0f - sx;
	float osy = 1.0f - sy;
	
	line(by)[bx] += c*osx * osy;
	line(by)[bxo] += c*sx * osy;
	line(byo)[bx] += c*osx * sy;
	line(byo)[bxo] += c*sx * sy;
}

void floatimage::AddPixelGaussian(float sx, float sy, float width, float alpha, fRGB c)
{

	int w	   = ceil(width);
	int xstart = max(0,sx-w);
	int xend   = min(sizex-1,sx+w);
	int ystart = max(0,sy-w);
	int yend   = min(sizey-1,sy+w);
	float expX = expf(-alpha * width * width);
	float expY = expf(-alpha * width * width);
	float total_weight = 0.f;

	for (int m = ystart; m < yend; m++)
		for (int n = xstart; n < xend; n++)
		{
			float xdist2 = (n-sx)*(n-sx);
			float ydist2 = (m-sy)*(m-sy);
			total_weight += max(0.f, float(expf(-alpha * xdist2)) - expX)*
				max(0.f, float(expf(-alpha * ydist2)) - expY);
		}

	for (int y = ystart; y < yend; y++)
		for (int x = xstart; x < xend; x++)
		{
			float xdist2 = (x-sx)*(x-sx);
			float ydist2 = (y-sy)*(y-sy);
			float weight = max(0.f, float(expf(-alpha * xdist2)) - expX)*
				max(0.f, float(expf(-alpha * ydist2)) - expY);

			line(y)[x] += c*(weight/total_weight);
		}
}

void floatimage::AddPixelGaussianNormalized(float sx, float sy, float width, float alpha, float N, fRGB c)
{

	int w	   = ceil(width);
	int xstart = max(0,sx-w);
	int xend   = min(sizex-1,sx+w);
	int ystart = max(0,sy-w);
	int yend   = min(sizey-1,sy+w);
	float expX = expf(-alpha * width * width);
	float expY = expf(-alpha * width * width);

	for (int y = ystart; y < yend; y++)
		for (int x = xstart; x < xend; x++)
		{
			float xdist2 = (x-sx)*(x-sx);
			float ydist2 = (y-sy)*(y-sy);
			float weight = max(0.f, float(expf(-alpha * xdist2)) - expX)*
				max(0.f, float(expf(-alpha * ydist2)) - expY);

			line(y)[x] += c*(weight/N);
		}
}

float floatimage::GetGaussianNormalized(float width, float alpha)
{

	int w	   = ceil(width);
	float expX = expf(-alpha * width * width);
	float expY = expf(-alpha * width * width);
	float total_weight = 0.f;

	for (int y = -w; y < w; y++)
		for (int x = -w; x < w; x++)
			total_weight += max(0.f, float(expf(-alpha * x*x)) - expX)*
				max(0.f, float(expf(-alpha * y*y)) - expY);

	return total_weight;
}

/*
void floatimage::Gaussian(float width, float alpha, fRGB c)
{
	int w	   = ceil(width);
	float expX = expf(-alpha * width * width);
	float expY = expf(-alpha * width * width);

	for (int y = 0; y < sizey; y++)
		for (int x = 0; x < sizex; x++)
			for(int i = max(0,y-w); i < min(sizey-1,y+w); i++)
				for(int j = max(0,x-w); j < min(sizex-1,x+w); j++)
				{
					float xdist = fabs(j);
					float ydist = fabs(i);
					line(y)[x] += c*max(0.f, float(expf(-alpha * xdist * xdist)) - expX)*
					max(0.f, float(expf(-alpha * ydist * ydist)) - expY);
		}
}
*/

void bubbleSort9(float x[]) ;

void floatimage::MedianFilter3x3()
{
	floatimage temp;
	temp.Allocate(sizex, sizey);
	
	float xr[9];
	float xg[9];
	float xb[9];
	int pos;
	
	for (int j = 0; j < sizey; j++) {
		for (int i = 0; i < sizex; i++) {
			if(i == 0 || j == 0 || i == sizex-1 || j == sizey-1)
				temp.line(j)[i] = fRGB(0.f,0.f,0.f);
			else
			{
				pos = 0;
				for(int y = j-1; y <= j+1; y++) {
					for(int x = i-1; x <= i+1; x++) {
						
						float rscale = 0.001914;
						float gscale = 0.001396;

						xr[pos] = line(y)[x].r - 1.5*rscale;
						xg[pos] = line(y)[x].g - 1.5*gscale;
						//xb[pos] = line(y)[x].b;
						pos++;
					}
				}
				bubbleSort9(xr);
				bubbleSort9(xg);
				float val = __max(0.0,xg[1]) + __max(0.0,xb[1]); 
				//bubbleSort9(xb);
				//temp.line(j)[i] = fRGB(xr[1],xg[1],xb[1]);
				temp.line(j)[i] = fRGB(val,val,val);
			}
		}
	}
	
	SwapData(temp);
}

void bubbleSort9(float x[]) 
{
	for (int j=1; j < 9; j++) { 
		for (int i=0; i < 9-j; i++) {
			if (x[i] > x[i+1]) {
				float temp = x[i]; 
				x[i] = x[i+1]; 
			   x[i+1] = temp;
           }
       }
   }
}

void floatimage::NormalizeToMax()
{
	float max = -1.0;

	for (int j = 0; j < sizey; j++) {
		for (int i = 0; i < sizex; i++) {
			if(line(j)[i].r > max)
			{
				max = line(j)[i].r;
			}
		}
	}

	for (int j = 0; j < sizey; j++) {
		for (int i = 0; i < sizex; i++) {
			line(j)[i] /= max;
		}
	}
}
#ifdef	FLOATIMAGE_USEWVT
int floatimage::LoadWVT(const char *filename, float image_gamma)
{
	unsigned char *tmp;
	int	originx, originy;
	int ok = WvtFileToRGB(filename, &tmp, &originx, &originy);

	int img_width = originx > originy ? originx : originy;
	int power = 1;
	while(power < img_width)
		power *= 2;
	img_width = power;

	InitGamma(image_gamma);

	Allocate(originx, originy);

	int	img_plane_size = img_width*img_width;
	unsigned char *cur_r = tmp;
	unsigned char *cur_g = cur_r + img_plane_size;
	unsigned char *cur_b = cur_g + img_plane_size;

	for (unsigned int y = 0; y < sizey; y++)
	{
		for (unsigned int x = 0; x < sizex; x++)
		{
			line(y)[x].Set(gamma_table[cur_r[x]], gamma_table[cur_g[x]], gamma_table[cur_b[x]]);
		}
		cur_r += img_width;
		cur_g += img_width;
		cur_b += img_width;
	}

/*
	Allocate(img_width, img_width);

	int	img_plane_size = img_width*img_width;
	unsigned char *cur_r = tmp;
	unsigned char *cur_g = cur_r + img_plane_size;
	unsigned char *cur_b = cur_g + img_plane_size;

	for (int y = 0; y < sizey; y++)
	{
		float *buf = &line(y)[0].r;
		for (int x = 0; x < sizex; x++)
		{
			buf[0] = cur_r[x]/100;
			buf[1] = cur_g[x]/100;
			buf[2] = cur_b[x]/100;

			buf += 3;
		}
		cur_r += img_width;
		cur_g += img_width;
		cur_b += img_width;
	}
*/
	return ok;
}

int floatimage::SaveWVT(const char *filename, float image_scale, float image_gamma)
{
	// buffer : buf
	// width  : sizex
	// height : sizey

	unsigned char lookup[65536];
	for (int x = 0; x < 65536; x++)
	{
		float tval = (float)(::pow((double) short2float(x) * image_scale,(double)  1.0 / image_gamma) * 256.f);
		lookup[x] = (tval > 255.0 ? 255 : (unsigned char) tval);
	}

	// get the nearest power-of-two square size
	int img_width = sizex > sizey ? sizex : sizey;
	int power = 1;
	while(power < img_width)
		power *= 2;
	img_width = power;
	int	img_plane_size = img_width*img_width;

	unsigned char *tmp = new unsigned char[img_plane_size * 3];
	unsigned char *cur_r = tmp;
	unsigned char *cur_g = cur_r + img_plane_size;
	unsigned char *cur_b = cur_g + img_plane_size;

	if (!tmp)
		return 0;
	else
		memset(tmp, 0, img_plane_size*3);

	// go through each pixel (in each channel)
	// deinterleave and color correction
	for (int y = 0; y < sizey; y++)
	{
		float *buf = &line(y)[0].r;
		for (x = 0; x < sizex; x++)
		{
			cur_r[x] = lookup[float2short(buf[0])];
			cur_g[x] = lookup[float2short(buf[1])];
			cur_b[x] = lookup[float2short(buf[2])];

			buf += 3;
		}
		cur_r += img_width;
		cur_g += img_width;
		cur_b += img_width;
	}
	
	// encoding
/*
	// debug use
	FILE *fp = fopen("test.raw", "wb");
	if(fp == NULL)
		return 0;
	fwrite(tmp, img_plane_size*3, 1, fp);
	fclose(fp);
*/
	int ok = RGBToWvtFile(filename, tmp, img_width, sizex, sizey);

	delete [] tmp;

	if (!ok) 
		return 0;
	else
		return 1;
}
#endif